CHEBYSHEV COLLOCATION FOR LINEAR, PERIODIC 
ORDINARY AND DELAY DIFFERENTIAL EQUATIONS: 
A POSTERIORI ESTIMATES 



ED BUELER 



Abstract. We present a Chebyshev collocation method for linear ODE and DDE prob- 
lems. Theorem in section [3] gives an a posteriori estimate for the accuracy of the 
approximate solution of a scalar ODE initial value problem. Examples of the success 
of the estimate are given. For linear, periodic DDEs with integer delays we define and 
discuss the monodromy operator U in section 01 Our main goal is reliable estimation of 
the stability of such DDEs. In sectional we prove theorem 1141 which gives a posteriori 
estimates for eigenvalues of U, our main result. Theorem 1141 is based on theorem 1111 
a generalization to operators on Hilbcrt spaces of the Bauer-Fike theorem for (ma- 
trix) eigenvalue perturbation problems. Section describes the generalization of these 
results to systems of DDEs and an example is given in section The computation of 
good bounds on ODE fundamental solutions is an important technical issue and an a 
posteriori method for such bounds is given. An additional technical issue is addressed 
in section |SJ namely the evaluation of polynomials and of the L°° norms of analytic 
functions. Generalization to the non- integer delays case is considered in sectional 
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1. Introduction 
Consider the following delay differential equation (DDE): 

(I) y(t) = ay(t) + {b + sin(37rf)) y(t - 2). 

Question: For what values of a, b G R is equation (pj stable in the sense 
that all solutions y(t) go to zero as t — > oo? 

For a, b e [—3,3] x [—2,4], the answer is given by figure ^ a stability chart. 
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Figure 1. A stability chart for equation with a (numerically- 
determined) stable point (a, b) = (—1.1, 1) indicated. 

Pictures like figure ^ of the stability of linear DDE are useful in applications. Examples 
include engineering design and control problems for systems described by DDEs 01211, 
but also applications in biology [21], for instance. 

For constant-coefficient cases the stability can be determined "by hand" using complex 
variable techniques but, in general, the stability of linear DDE systems must be deter- 
mined approximately. If the Floquet transition matrix for a linear, periodic ordinary 
differential equation (ODE) associated to a linear, periodic DDE with integer delays can 
be found exactly, then complex variable techniques will determine stability of the DDE 
(see section 8.3 of JI]). In general, however, ODE Floquet matrices must themselves 
be approximated. The techniques of this paper directly approximate DDE stability by 
directly approximating the DDE monodromy operator (section HJ) associated to equation 



CHEBYSHEV COLLOCATION ON ODES AND DDES 



3 



Figure Q was produced pixel-by-pixel by approximating the monodromy operator for 
the parameter values a, b associated to each pixel in the chart, and (numerically) com- 
puting the largest eigenvalue of the approximating matrix, a monodromy matrix. The 
approximation method can be described as Chebyshev (spectral) collocation. 

The first question about this procedure, mathematically, is how accurately one can 
solve initial value problems for DDE (^)? That is, how accurately can one perform the 
"method-of-steps" JT] for the DDE? This paper proves a posteriori estimates for the 
accuracy of the collocation algorithm on such problems. 

The next question, mathematically, is about the quality of the eigenvalues of the mon- 
odromy matrix as estimates of the eigenvalues of the monodromy operator. The mon- 
odromy operator is compact and thus its spectrum is entirely eigenvalue spectrum (appen- 
dix B). In this paper we combine the results of section El with an eigenvalue perturbation 
theorem for operators on Hilbert spaces (theorem ^2 in section EJ to give a posteriori 
estimates for the distance between a nonzero eigenvalue of the monodromy operator and 
the nearest eigenvalue of the monodromy matrix. 

For instance, suppose a = —1.1 and b = 1 in equation (JTJ), a point indicated in figureHJ 
to be stable but near the stability boundary. Figure El shows as solid dots the computed 
eigenvalues of the monodromy matrix. For the largest three eigenvalues, discs are shown 
whose radii r = .0434 are proven maximum distances by which the large eigenvalues of 
the monodromy operator might differ; see theorem El By "large eigenvalues" we mean 
that we have chosen to consider eigenvalues /jgC such that \fi\ > 5 = 0.2. 

Because the largest eigenvalue has magnitude = .9369 and thus + r < 1, 
figure El represents a proof that the parameters (a, b) = (—1.1, 1) are eigenvalue stable for 
equation (JTJ). 

The numerical eigenvalues in figureEl the dots, were generated using a (iV+ 1) x (iV+ 1) 
monodromy matrix with N = 184. It turns out that the radius of the discs decays 
roughly exponentially with N as illustrated in figure 01 In this example, when N ^ 220 
we find that floating-point errors, and the fact that the condition number of the matrix 
of eigenvectors of the monodromy matrix is roughly 10 6 , limits accuracy to roughly 10 -5 . 

Subsection 15. 5l presents this scalar DDE example in detail. Section[7|addresses a system 
of DDEs, a delayed Mathieu equation, with similar results. 

This paper is somewhat more expository than is standard in the numerical analysis 
literature. We start by recalling the basics of Chebyshev (spectral) collocation in section 
El Then we prove new results for the first order scalar DDE case in sections 0] and 
El We then generalize the same machinery to systems of DDEs (section |HJ) and apply 
the machinery to the example of a delayed, damped Mathieu equation (sectional). (In 
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Figure 2. If a = —1.1 and 6=1 then the largest three eigenvalues of the 
monodromy operator for equation (P) are proven to be within the given 
discs. Thus the parameter point (a, b) is stable. (Unit circle and circle 
r = 5 = 0.2 also shown.) 
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FIGURE 3. For DDE ((H) with (a, b) = (-1.1,1): The monodromy matrix 
is an (N + 1) x (N + 1) matrix approximation to the oo x oo monodromy 
operator. The radius found from theorem El of the error discs around the 
numerical eigenvalues of the monodromy matrix decays exponentially with 
N. 



a less expository paper we would do the systems case immediately, with inevitable loss 
of clarity to the non-expert.) Further generalizations and technical issues are covered in 
sections |H1 and El and in the appendices. 
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We now review some of the numerical analysis literature relating to stability of DDEs 
by spectral methods. Spectral methods were perhaps first applied to DDE initial value 
problems in |3 j and though in the former work the emphasis is on "/i-refinement" 
and in the latter the techniques are limited to constant-coefficient linear problems. More 
recently, collocation using the Gauss- Legendre points, in particular, is applied to the 
problem of finding periodic solutions to nonlinear DDEs in [5] and the stability 1 and 
convergence under /i-refinement of the resulting method is addressed in jHj. 

On the other hand, there exists some engineering literature addressing stability (as 
opposed to solutions of initial value problems) of ODEs and DDEs by spectral methods. 
For ODE problems we specifically note [TH] who use equally-spaced collocation points and 
[23 [23] who use a Chebyshev Galerkin (or Galerkin-collocation, respectively) method. In 
fact, such engineering applications led the current author and coworkers to DDE stability 
questions [Hj. 

The current author knows of two numerical analysis papers, namely JH] and [201, which 
address the application of spectral methods to the stability of linear DDEs and functional 
differential equations. The former surveys collocation methods for approximating eigen- 
values of monodromy operators (i.e. Floquet multipliers) and gives numerical evidence 
for spectral convergence (see table 2 of JH])- I n the latter work it is proven that the 
numerical method preserves a property ("RHP-stability") which suffices for stability of 
the functional equation. In neither case are estimates on the eigenvalues addressed. 

This paper grows out of enjoyable work on DDEs in engineering applications initi- 
ated and sustained by Eric Butcher, and continued with students Victoria Averina, Tim 
Carlson, Haitao Ma, Praveen Nindujarla, Jake Stroh, and Ben White. 

2. Interpolation and spectral differentiation at Chebyshev points 

We are interested in solving ODEs or DDEs on a fixed interval t E I and we suppose 
I — [—1,1] as this choice is convenient relative to Chebyshev polynomial conventions. All 
results are easily translatable to other intervals. The functions we consider will generally 
be continuous and complex- valued; denote the space of such functions by C(I). (It will 
be most useful to consider complex-valued functions when considering stability problems 
as the eigenf unctions of the real monodromy operator (section [3]) are generally complex.) 
We use the norm ||/||oo = max ie / |/(^)| for / G C(I). 

Definitions. Let N > 1. Let Vn be the space of at most degree N polynomials with 
complex coefficients. Let Cn = {to, ■ ■ ■ , ^v} C / be the Chebyshev collocation points 
(extreme points) tj = cos(nj /N). Note t = 1, t N = —1, and < tj. 

1 We distinguish here between stability of the DDE and stability of the numerical method. 
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As a subspace of C(I) the dimension of Vn is N + 1. For instance, the monomials 
l,t, t 2 , . . . ,t N are a basis of Vn, as are the Lagrange polynomials Lj(t) = Ylk=o k^jft ~ 
tk)/(tj — tk) for the collocation points. 

Polynomial interpolation using the collocation points is our first concern. (See [3] or 
[22] for the basic theory of polynomial interpolation.) Consider the following maps: 

Definitions. 

• Let £n '■ C(I) — > C N+1 be the linear map which evaluates the input function at 
the m collocation points tj, so (£jv/)j = f(tj) for j = 0, . . . , N. 

• Let E N : Vn — > C N+1 be the restriction of £n to polynomials. 

• Let P N : C N+1 — > Vn be the map which gives the unique polynomial of degree 
at most N = m — 1 having the input values at the collocation points tj. Note 
(Pjyy)(t) = J2f=o v jLj(t), though we will have no practical need for this formula. 

• Let In = Pn&n '■ C(I) Vn be the interpolation map which constructs the 
unique polynomial of degree iV which agrees with the input function at the collo- 
cation points. 

Note that E N , Pn are linear isomorphisms, P^ 1 = E N , and In is a projection (i.e. I N = 
In). 

The quality of polynomial interpolation depends on the regularity of the function being 
interpolated. Nonetheless it is possible to find analytic functions f on I and (families of) 
interpolation points such that / is not well-interpolated by polynomials • Fortunately, 
Chebyshev collocation points are excellent for interpolation in the analytic case: 

Theorem 1 (proven in |2Z|, in particular). Suppose f is analytic in an open, simply- 
connected region R such that [—1, 1] C R C C Then there exist constants C > and 
p > 1 independent of N such that if p — In} £ Vn is the polynomial interpolant of f at 
the N + 1 Chebyshev collocation points Cn, 

(2) Wf-vWoo^Cp^. 

Furthermore, suppose E C R is an ellipse with foci ±1 — see figure^ below. If S, s are the 
lengths of the semimajor, semiminor axes of E, respectively, then inequality (J2J applies 
with p = S + s. In fact, for k = 0, 1, 2, . . . there exists Ck > such that 

\\f {k) -p {k) \\oo<C k p- N : 

for p = S + s, where f^ is the kth derivative of f on I. 

From polynomial interpolation one may construct spectral differentiation. Let 

D N = En^-Pn : £ N+l - C^ + \ 
at 
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Figure 4. If / is analytic in a region EcC which contains a "regularity 
ellipse" E with foci ±1 and semiaxes S, s then Chebyshev interpolation 
converges at exponential rate p~ x where p = S + s. Note that S 2 = 1 + s 2 
for an ellipse with foci ±1. 

whose (N + 1) x (N + 1) matrix with respect to the standard basis for C N+1 is called 
the Chebyshev spectral differentiation matrix. The matrix "differentiates" v G C N+1 
as follows: construct a polynomial p with values v at the collocation points; differentiate 
it p' = and evaluate the result at the collocation points (D N v)j = p'(tj). See (2H] for 
explicit formulas for the entries of and applications of D^. Note that all entries of 
.Djv are real and that is nilpotent. 

Exponential decay of error as iV increases is known as spectral convergence, so theorem 
[U says polynomial interpolation at the Chebyshev collocation points yields spectral con- 
vergence for analytic functions. The goal in this note is to show spectral convergence for 
the solutions of ODEs and DDEs, and for certain eigenvalue problems, by the numerical 
methods described below. 

We are interested in linear ODEs and DDEs which have highly regular (or piecewise- 
regular) coefficient functions. In particular, because this assumption applies to the engi- 
neering applications which we have addressed and/or foresee jH], we suppose these coef- 
ficient functions are (piecewise) analytic. In the DDE case we will additionally suppose 
that history functions are (piecewise) analytic. 

We end this section with a technical tool, namely, polynomial interpolation by discrete 
Chebyshev series. Recall that Tf-(t) = cos(A;arccost) are the Chebyshev polynomials 
and that the collocation points tj = cos(jir/N) are the extreme points of T^(t) since 
T N (tj) = cos(Njn/N) = . If f E C(I) then p = I N (f) is an N degree polynomial 
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which may be found by the following formulas 

(3) p(t) = ^/ fc T fe (t), fk^Cvfte), C kj = ——cos(^-), 

k=0 j=0 ^ j ^ k ^ ' 

where 7^ = 2 if j = 0, iV and 7j = 1 otherwise. These formulas may be implemented 
by a modification of the fast Fourier transform |29j . That is, the FFT may be used to 
compute P/v and £jv if desired. 

3. A posteriori ESTIMATES FOR INITIAL VALUE PROBLEMS 

3.1. An algorithm and error estimate. Consider the scalar, linear ODE initial value 
problem 

(4) y(t) = a(t)y(t)+u(t), y(-l) = y . 

Needless to say, this problem has a well-known exact solution. Let $ a (£) be the ODE 
fundamental solution solving $ a (/) = a,(t)<& a (t) and $ a (— 1) = 1. In this scalar case 
&a(t) = exp ^ a(s) ds^j, of course, but the existence of the fundamental solution $ a (£) 
generalizes to the systems case while this exponential formula does not (at least, directly 
). The exact solution to (HJ) is then 



(5) y{t) = 



Vo + J 1 u(s)ds 



by variation-of-parameters. 

Let Z)jy be the (N + 1) x (N + 1) matrix modification of which satisfies (-Djv)^ = 
(D^)ij for rows i — 1, . . . , N and all columns, but has last row 

pjv)jv+i,fe = 0, 1 < k < N, and (£>n)n+i,n+i = 1- 

The product D N v represents both "?/(£)" for — 1 < t < 1 and u y(— 1)" if v is the (Cheby- 
shev) collocation approximation of y{t). For a 6 C(I) let M a be the (AT + 1) x (N + 1) 
diagonal matrix with entries 

(M a )i, 

The product M a v represents the multiplication u a(t)y(t)" for —1 < t < 1 if v is the 
collocation approximation of y{t). For u G C(I) and ?/o £ C let u G C Ar+1 have entries 

(6) Ui = 
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That is, Ui = {£^u)i if 1 < i < N while un+i = Ho- The following is immediate from the 
construction of L>n, M a , and u. 

Lemma 2 (the collocation algorithm). Let a, u G C(I) and let y G C T/ie following 
two descriptions of the collocation approximation of the problem (0J) are equivalent: 

• p G "Pat satisfies 

(7) p(tj) = a(*iM*j)+«(*i). 0<j<iV-l, and p(-l) = y ; 

• v G C^ 1 satisfies 

(8) D N v = M a v + u. 
The equivalence is p = Pjyv and v = E N p. 

We will call the numerical algorithm described by 

y(t) w p(t) = = P^(,D W - Ma)" 1 ?}, 

that is, by the above lemma, the collocation algorithm. (For a complete numerical algo- 
rithm we would need to specify techniques for constructing and for solving the linear 
equations which arise, but for now we assume these tasks are done exactly. See, however, 
section |H1) 

The algorithm generalizes gracefully to scalar DDEs (sections |U and EJ) and to system 
of ODEs and DDEs. It is well-proven in practice. 

Our goal of showing spectral convergence for the collocation algorithm is, roughly 
speaking, the goal of showing that the collocation algorithm does as good a job of ap- 
proximating the solution to the initial value problem as would polynomial interpolation 
on y{t) in ©. 

What we can prove is an estimate on the maximum difference between y(t) and p(t) in 
terms of computable quantities. 

Theorem 3. Suppose y G C(I) solves initial value problem (|4j) with a G C(I). Suppose 
p G Vn, N > 1, is calculated by the collocation algorithm (lemma&) . Let 

C a = exp max{Re a(s), 0} ds^j 

and 

R p = p(-l) - a(-l)y - u(-l). 
Then the uniform error of p as an approximation to y satisfies 

(9) \\y -p||oo < 2E a (\\ap - InMWoo + \\u - /at(w)||oo + |-Rp|)- 
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Note that |$ (t)| = exp ( Rea(s) dsj < C a , so C a simply bounds the fundamental 
solution. The quantity Rp is a measure of the degree to which the approximating poly- 
nomial p does not satisfy the differential equation at the initial time; it could be called 
the initial residual for the method. Actually computing the initial residual is easy, as 
p(-l) = (D N v) N+1 . 

It would be nice to show that \\ap — Jjv(ap)||oo an d R P are (exponentially) small in 
iV and thus make this theorem into a (spectral) convergence result for the collocation 
algorithm. For now, however, this theorem is "a posteriori" because one must actually 
compute p in order to evaluate the error estimate. 

Computing the maximum norms \\ap — I^{ap) ||oo, \\u — /jv(^)||oo is done by the tech- 
niques described in section |8] 

In the special case where a(t) = Oq is a constant, we can be more specific about the 
conclusion of this theorem. If Re do > then 

(10) \\y-p\\oo <2e 2Rcao \\u-I N {u)\\ oc + c\R p \ 

where c= (vr(|a | + 1) + 4) e 2Rea °/(2 N 2 ) and R p = p(-l) - a y ~ u(-l). If, Rea < 0, 
on the other hand, then 

(11) \\y-p\\oo < 2||tt-Ijv(u)|| (XJ + min|c,-^} \R p \, 

here c = (vr(|a | + 1) + 4) e~ 2Rea °/(2 iV 2 ). Note that the constant multiplying the initial 
residual actually decays with increasing N. These inequalities are proven in appendix A. 

The error estimate in theorem El is absolute, not relative. In practice one may choose 
iV so that \\y — p||oo <C ||p||oo, established using the estimate, and then increase N, if 
needed, so that \\y — p\\oo/\\p\\oo < 5 for some relative tolerance 5. In the case that \\p\\oo 
is itself small, absolute error \\y — pW^ is exactly what we would wish to control. 

3.2. Examples. Before proving this theorem we wish to demonstrate its effectiveness in 
simple examples. In particular, the behavior of different cases a(t) = ao (exponential 
growth when Reao 3> 0, decay when Rea <C 0, and combined with oscillation when 
| Ima | is large) is captured by equations (fTUj) and (JTTJ). We start with an easy case. 

Example 1. Consider the following initial value problem 

y{t) = 3y{t)+t, y(-l) = y , 

with exact solution y(t) = e 3 (' +1 ) (y — |) — | (t + |). We choose y = 0.2 so that y{t) 
takes values in the interval [—9.5, 0.2] for t G J and thus y(t) is 0(1) on I. Since u(t) = t 
is a polynomial of degree one and a(t) = 3 is constant, special case (fTU|) gives the estimate 

\\y - Plloo < 2(tt + l)e 6 N~ 2 \p(-l) + 0.4| 
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for the maximum difference between the collocation polynomial p e Vn and the exact 
solution y. We compare this estimate to the actual maximum error, computed by sampling 
the exact solution and approximating polynomial at 1000 equally-spaced points in I. See 
figure El Exponential decay of the actual error is clear. The convergence stops roughly 
3 orders of magnitude above machine precision e m ~ 2 x 10 -16 . Note that \\y\\oo ~ 9.4 
so the relative error \\y — p\\oo/\\y\\oo is roughly two orders of magnitude above e m . The 
estimate, however, is consistently only about about a factor of 10 too large, and this is 
quite acceptable because the convergence is so fast that the estimated value for iV to give 
a certain error is too big by only one in practice. 2 
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Figure 5. Actual error — p||oo and the estimate (fTUJ) for the initial value 
problem in example ^ 



Example 2. Cases with rapid exponential growth, decay, or oscillation — i.e. stiffness — 
show that spectral methods have to work hard, too. Let u(t) = and consider y(t) = 
ay(t), y(—l) = 1, with exact solution y(t) = e a( - t+1 \ Figure |B1 shows actual error and the 
estimate in cases a = 10 and a = —10. The behavior of the actual error is quite different 
in the two cases, and we use special cases (fTD|) and (fTTj) . respectively. For a = 10, note 
that the solution y(t) = e 10( * +1 ) grows from y(—l) = 1 to y(l) = e 20 ~ 5 x 10 8 , and 
thus the error is reasonable in a relative sense. Also, convergence does not start until 
iV = 20 or so. The significance of this number is that a polynomial of degree at least 20 
is required to approximate the rapid growth of e 10t on I. For a = —10 we see different 

2 The estimate might be suspiciously good — it tracks the behavior of the actual error when both are 
comparable to machine precision — until one recalls the a posteriori nature of the estimate. 
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behavior, with immediate but slower convergence. The explanation is that low order 
polynomials do exist which approximate e~ 10 ' in an absolute sense, but only polynomials 
of high degree succeed in being comparably small to e~ 10i for t G / close to +1. For both 
a = ±10 spectral convergence ends at about iV = 30, as non-truncation numerical errors 
(e.g. errors in solving linear equations and rounding errors) become dominant. 
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Figure 6. (a) Estimated and actual error for y = lOy. (b) For y = —lOy. 

Example 3. The nonhomogeneous term u(t) is usually not a polynomial. Suppose a(t) = 
a Q = 3 + 37i, y(—l) = 0.2, and u(t) = sin(20£). The exact solution is found by variation- 
of-parameters. We calculate the ||w — In(u)\\oo term in theorem El by evaluating u and 
In{u) at 1000 equally spaced points. Figure [7| shows four quantities: the estimated error, 
the actual error, the initial time residual error \p(— 1) — ayo— u{— 1)|, and the interpolation 
error ||u — Iat(w)||oo- The need for N to be sufficiently large to "handle" the oscillations 
present in $ a (t) = e~ ao( - t+1 ^ and u{t) is clear. 

Example 4. Finally, the leading coefficient a(t) need not be constant. Suppose a(t) = 2t, 
u(t) = tsin(3t 2 ), y(— 1) = 1. The exact solution can be found by variation of pa- 
rameters because $ a (£) = e* 2_1 so that the relevant integral is f e l ~ s2 s sin(3s 2 ) ds = 
e J e~ z sin(3z) dz (under the substitution z = s 2 ). Note C a = e is a conservative bound 
on $ a (t). Figure |H1 shows that though the initial residual error happens to be small for all 
N, the interpolation errors \\ap — ijv( a P)||ooj ll M — In(u)\\oc are both active contributors to 
the error estimate and the actual error \\y — p\\oo- These quantities decay exponentially 
till error and estimate are of close to machine precision 10~ 15 . 

The collocation algorithm works well on initial value problems with or without the 
estimate given in theorem 01 including on ODE (and DDEs) systems and in cases with 
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Figure 8. Estimated and actual error for y = 2t y + 1 sin(3t 2 ) , y{— 1) = 1. 
Both interpolation error terms on the right side of inequality (JHJ) contribute 
to the estimate. 



multiple delays and merely piecewise analytic coefficients. We generally extract roughly 
double precision results — say, 8 to 14 correct figures for ||a||oo ^ 10 — in those problems 
for which we can compare to exact answers. 
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A caveat: Theorem El assumes that the spectral differentiation matrix Dn is exactly 
computed and that linear system (JBJ) is exactly solved. Computation with large N (with 
iV « 100 already fairly large) must be carefully handled to extend spectral convergence 
or even to maintain stability (see section |H] and references and [3*U]\ 

3.3. Proofs. We now turn to proving theorem |3J First we establish a lemma which is 
reminiscent of computations with the Dirichlet kernel in Fourier series |16j . namely, we 
rewrite the monic polynomial /jv(t) with roots to, ... , tjv-i as a closed-form trigonometric 
expression. Figure El shows graphs of the polynomials Zjv(t) f° r some values of N. 



Lemma 4. For N > 1 let Zjv(t) = (t — t )...(t — tjv-i)- -For t e I and t = cosO, 

(cos(0) - 1) sm(N 9) 
sin(0) 

and Zjv(-I) = (-l) iV A r 2 2 - iV ; IMU = N 2 2 ~ N , in particular. 



WW = ojv-i^/m ' -Kt<l, 




-1 -0.5 0.5 1 



t 

Figure 9. Polynomials Zjv(^) for three values of iV. It turns out that Zjv(^) 
has a closed-form expression in = arccost. The dots at t = — 1 have 
values (-l) 7V iV2 2 - Jv . 



Proof. We use discrete Chebyshev series (section |2J) to write /^(t) as a sum of Chebyshev 
polynomials T k (t). Noting that /^(i,) = for j = 0, . . . , iV — 1, 

(12) l N (t) = [1 - 2T.it) + 2T 2 (t) - • ■ ■ + (-l^-^T^t) + (-l) N T N (t)} . 

Now, since Jjv(t) is monic and because the coefficient of £ in T^(t) is 2 Ar_1 (use 
the recursion relation J* = 2tTfc_i — Tk-2), it follows that 1 = ^fe^ (— 1) JV 2 JV ~ 1 or 
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In(— 1) = (— l) N N 2 2 ~ N as claimed. From equation (fT2"|) we have 
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2 iV-l 



1 + 2^(-l) fc cos(£;fl) + (-l) Ar cos(^) 



fe=i 



We produce a closed form expression for geometric sum, 

N-l N-l 

l + 2^(-l) fc cos(^)= (- el6 ) k = (- r 

k=l k=-(N-l) 

and find, through trigonometric identities, that 



_ 1 cos((A^- l/2)g) 
cos(fl/2) ' 



In if) 



2 N-1 
1 

oJV-i 



cos(6y2) 
-sin(JV0) sin(0/2)' 
cos(fl/2) 



(cos(0) - 1) sin(iV0) 



)JV-1 



sin(0) 



□ 



Proof of theorem^ Let g = I^{ap) and w = Jatu. Note g(t,) = a(tj)p(tj) and u>(i/) = 
for j = 0, . . . , AT, of course. 
Let r = p — q — w. Note r £ and r(ij) = for j — 0, . . . , N — 1 by lemma 121 Thus 
there is /3 G C such that r{t) = (3l N (t). By lemmaEJand evaluating r(t) with £ = —1, 



(13) 



(3l N {~l) = p- 



1) N N 



)JV-2 



p(-l)-a(-l)p(-l)-u(-l) = i2p. 



On the other hand, noting p = q + w + (3 In = ap + (q — ap + w + (3 In) , variation-of 
parameters gives 

(14) V {t) = $ a (t) f 
Subtracting equations © and (jTljl yields 

(15) y {t)-p{t)= [ ¥ 8 [a{s)p{s) -q{s)+u(s)-w(s)]ds- /3 I ^ t s l N {s)ds. 



Vo + / $a0) 1 - a(s)p(s) + iu(s) + /3 l N (s)) ds 



where = $ a (t)$ a (s)" 1 . Note that |$*| < C a for s < t. Thus by lemmalU 

( • N 

\y(t) -p(t)\ < 2C a \\ap - q\\oo + \\u - w\\oo + 



)N-2 



Taking absolute values in equation (fT3*|) we have proven the theorem. 



□ 



Finally, a technical detail remains, needed when estimating eigenvalues of a DDE mon- 
odromy operator. 
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Corollary 5. Under the hypotheses of theorem^ we have the derivative estimate 

\\y-p\\°c < (2C a ||a|| 00 + l)^||ap-/j V (ap)|| 00 + ||« - /jv(«)IU + |-Rp|)- 
Proof. Since y = ay + u, p = q + w + (31 n, using the notation of the above proof we have 

y - p = a(y - p) + (ap - q) + (u - w) + (3l N 

so 

\\V -p\\oo < N|oo||Z/ ~P\\oo + \\aP - Q\\oc + \\U - HL + |-Rp|- 

The result follows from theorem 01 □ 
4. The monodromy operator for linear, periodic DDEs with integer 

DELAYS 

4.1. The monodromy operator U on C(I). Consider the scalar linear DDE problem 
for t G I: 

(16) x(t) = a{t)x{t) + b(t)x{t - 2), x(s) = f(s + 2) for - 3 < s < -1. 

We suppose a,b, f G C(J) and we extend a, 6 periodically to R. It is easy to solve problem 
(fT7)J) by variation-of-parameters to find x(t) corresponding to a given history function f(t). 
That is, for t £ I we simply solve the ODE initial value problem x(t) = a(t)x(t) + u(t), 
x(—l) = /(l) with w(t) = b(t)f(t). Suppose, furthermore, that a(t), b(t) are periodic with 
period T = 2. The variation-of-parameters formula for (|16J) can be translated to each 
interval of length 2 to find x(t) for all t > — 1, and we regard the effect as an operator on 
C(I). 

Definition. Define the DDE monodromy operator 3 U : C(I) — > C(J) by 



(17) ([//)(*) = $ B (t) 



/(1)+ / ^.(aJ-^W/Ca) ds 



-i 



In appendix B we show that U is a compact operator and thus, in some sense, is as 
much like a matrix as any infinite rank operator. 

The solution x(t) can be found "by steps" with x n G C(I), n > 0: 

X n+ \ = Ux nj Xq = f. 

The actual solution x(t), t > —1, is built by "concatenating" the steps: x(t) = x n (t — 
2n-2) if t G [2n-3,2n-l]. 

DDE dHH) is therefor stable if and only if U k f -> as fc -> c» for all / G C(J). 
Equivalently, if p(t/) = max Agcr (fj) |A| is the spectral radius of U [IT], p(U) < 1 if and 

3 Or perhaps the delayed Floquet transition operator. 
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only if U k f — > as k — > oo for all / G C(I). This is our stability criterion, an eigenvalue 
criterion though there exist more sophisticated measures |28j . 

We are thus motivated to consider eigenfunctions of U. Suppose Ux = Xx for A ^ 0. 
Then y = Xx solves the ODE y = ay + bx, that is, 

'18) x = ( a + —b] x. 



X 

In this scalar case the solution is 

(19) x(t) = x(-l) exp (J a(s) + jb(s) ds 

Of course, such a formula does not help us find A or i immediately, but we will use it to 
address the quality of polynomial interpolation of x in the next section. Note, however, 
that since (Uf)(— 1) = f(l), X = exp(a + b/X) where a = a(s) ds and b = b(s) ds. 
This allows determination of eigenvalues by complex variable methods. Such methods do 
not generalize to systems, however, because the fundamental solution $ a ,b,\(t) to equation 
()18j) is not generally found by straightforward integration [14 . 

4.2. A useful Hilbert space H . Theorem 1111 below, which will allow us to compare 
the eigenvalues of U to those of Un, requires operators acting on a Hilbert space. The 
obvious space is the one for which the Chebyshev polynomials form an orthogonal basis. 

Definition. Let L 2 = L^(I) be the Hilbert space of complex-valued (measurable) func- 
tions / on I such that j\ \f{t)\ 2 {\ - t 2 )~ 1/2 dt < oo. (Where "T" is for "Chebyshev".) 
Let 



(f,g)» = I f{t)g(t) (i-ty 1/2 dt, \\f\\l = (fj) 



L 2 



Definition. Define the (normalized) Chebyshev polynomials T (t) = (1 / ^/7r)T (t) = 

arccost). The set {Tk}^L is an orthonormal 
("ON") basis of L 2 . For / e L 2 let fk = (Tk, jJ ^ the Chebyshev expansion coefficients 

of /, so fit) = J2T=ofkTk(t) (with convergence in L 2 ) and ||/|| 2 2 = Y,T=o \fk\ 2 - 

We might hope that U is a nice operator on L 2 with matrix entries we can calculate. 
Unfortunately, U is not even bounded on L 2 . This is because formula (|17j) defining U 
refers to the pointwise values of /, while elements of L 2 are actually equivalence classes 
which can take any values on given sets of measure zero. 

Exercise 1. Assume b(t) = 0. Construct a sequence of functions /„ G C(I) such that 
/ n (l) = 1 and \\f n \\ L 2 -> but \\Uf n \\ L 2 0. 
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Thus we need to refine our choice of Hilbert space. 
Definition. We define H l = H^(I) in terms of the expansion coefficients, 

H 1 = If e L 2 |f> + *0Ul 2 <ool c L\ 

I fc=0 J 

For /, g G H 1 we define the inner product and norm 



(f,g) Hl = £(1 + k) 2 f k g k , \\f\\ 2 Hl = (f,f) Hl . 



k=0 

Let 



f k (t) = (l + k)~ 1 f k (t) 



y/lft T (t), k = 0, 

v ^A(l + A;)- 1 r fc (t), k>l, 



the if ^normalized Chebyshev polynomials. The set {Tfc}^ is an ON basis of H 1 . 

The Hilbert space H l is introduced in [21] . It is a Sobolev space but it is not equivalent 
to W\ = |/ G L 2 1 1|/ 11^2 + ||/||l2 < oo|. Elements of H 1 are, more regular that those of 
L 2 . In fact, H 1 C C(I) and pointwise evaluation is bounded. 

Lemma 6. If f G i/ 1 i/ien / is continuous (has a continuous representative f ) and 
\f(t)\ < 0.9062||/||^i. In particular, 

oo oo „i 

(20) 5J = fkT k (l) = E / ^(^X 1 " * 

fe=0 fc=0 



zs a bounded linear functional on f with \\5i\\ < 0.9062 and 8\f = /(l). 

Proof. The sum defining t^Z converges absolutely, and, in fact, Y^=a fkT k {t) converges 
uniformly for f G H l : 



OO OO OO 

v^E lAifcwi < E = E YTi ' (1 + 

k=l k=l k=l 

/ oo \ V2 / °° \ V2 

< E( 1+ *) _a E^+^W - 



CHEBYSHEV COLLOCATION ON ODES AND DDES 19 

using Cauchy-Schwarz. It follows that / is continuous because (in addition) = 
/ a.e. Note ££° =1 (1 + k)~ 2 = tt 2 /6 - 1. Thus 



fc=0 \fe=l 

1/2 



1 V2 

|2 



< 



(2/tt)|/ | 2 + 2C 2 £(1 + fc) 2 | f k \ 2 < V2C\\f\\ m 



k=i 



where C 2 = vr/3 - 2/vr and using a + o < V2a 2 + 2b 2 . Note v^C < 0.9062. □ 

Equation flUJ) is meant to suggest "<5i(t) = £^T*(1)T*(*)(1 - t 2 ) _1/2 " as a "delta 
function" , but this sum does not converge. 

In using the Hilbert space H 1 we will need to know that if / is sufficiently regular on 
/ then / G H . We give a criteria via the Fourier series of f(cos8). 

Definition. Let L 2 F be the space of measurable functions g on [— ir, n] such that J™ \g(0) | 2 d6 < 
oo. (Here "F" is for "Fourier".) Let \\g\\ 2 F = f^\g(6)\ 2 d6. 

Lemma 7. Suppose f G C 1 (J). The even function f{6) = /(cos#) is in Cp er [— 7T, ir], 
that is, f can be periodically extended with period 2it to be C 1 on IR. For k G Z let 
}(k) be the kth Fourier coefficient of f, that is, f_(k) = (2n)-^ 2 e- lke f{6) d6. Then 
fk = Hk) = K-k) fork>0 and f = 2-^1(0) . 
It follows that 

1 00 1 
(21) 11/11^ = 2 E l/(*)l a =2ll/l^ and 

k=—oo 



(22) ||/f ffl = i Ea + *0 2 l/(*)l a < II/IIf + ll/'lll- 

fc=— oo 

Thus C\I) CH 1 C C(J) and 

(23) H/Hk <2 7 r||/|| 2 + 2vr||/|| 2 

iffeC\i). 

Proof. The first assertion is elementary calculus. Now 

fk = ( f k,f) L2 = ^J\os(ke)i(e) de = \\J\f_ cosmm de 
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if k > because / is even. Because / is even it also follows that f(k) is real and thus 
fk — f(k). Similarly for f . Using Parseval's equality, equation (j21j) follows from 



E \l( k )\ 2 = ii(°)i 2 + 2 E toi 2 = 2 E lAi 

fc=— oo fc=l fc=0 



Next, if # G C£ er [-7r, tt] then ||^||| = JX-oc ^Wl 2 < 00 > and thus ©• Finally, 
since f'{9) = — f (cos 9) sin 9 , inequality (J22J) gives 

ll/ll^< r \f(cos9)\ 2 d9+ \f(cos9)\\m 2 9 d9<2n\\f 111 + 2*11/111. 

J— 7T J —7T 

□ 

4.3. The monodromy operator £7 on if 1 . We can now define the monodromy op- 
erator acting on H l . We then find its norm and estimate the accuracy of polynomial 
interpolation of its eigenfunctions, crucial steps in using theorem below. 

Definition. Suppose a,b e C(I) in DDE (JH|. Define U G C(H l ) by 
(24) ((if) (t) = $ a (t) f(l) + l\(s)- 1 b(s)f(s)ds 

for / G H 1 . 

The eigenvectors of £7 and £7 are identical. In fact, let in\ : H 1 C(I) be the injection 
and note: 

Lemma 8. £7, U are related by U = Uoi H i on H 1 and furthermore Uf = Xf for f G C(I) 
and X ^ if and only if f G H 1 and Uf = Xf . 

Proof. Note that if / G C(I) then £7/ G C 1 ^) C H 1 by definition of £7 (formula (HU)) and 
lemma But if £7/ = Xf it follows that / G H 1 and thus £7/ = Uf = Xf. Conversely, 
if / G H 1 then / G C(f) by lemma El and thus £7/ = Uf. □ 

In appendix B we see that the compactness of U follows from that of £7. On the other 
hand we have an a priori estimate of the norm of £7. 



Lemma 9. If C a = exp ( max{Re a(s), 0} ds 



then 



2 



as an operator on H 1 , where Cq = ||£>||oo; c i = 2.3(1 + ||a||oo) + 7r||&||oo; an d c 2 

7T>/2||a||oo||&||oo- 
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Proof. Clearly \\Uf\\ m < \\<S> a \\m\5if\ + \\$ a g\\m where g(t) = j t _ l <t> a (s)- l b(s)f(s)ds. 
Note |$ (t)| < C a and |6 a (t)| < C |H|oo- (Also ^(^(s)- 1 ! < C a if s < t; this is 



useful in below.) Thus \\<$> a \\m < V^C a (1 + 



12 \l/2 



by (|2SJ) - It follows that 



l^all^l^/l < 2.2715(1 + ||o|| 0O )C a ||/|| H i 

by lemma |H1 and using (a 2 + 5 2 ) 1 / 2 < a + 5 for a, 6 > 0. 
Now note 

(25) max{\^ a (t)g(t)\,\g(t)\}<C a \\b\\ c 



\f(s)\ds< VV2C a | 



L 2 



-1 



using Cauchy-Schwarz. Again from (J22J), noting $ a (t)g(t) = b(t)f(t) 

\\Qag\\h < 2vr||$ a ^||L + 27T||<i> a £ + $ a ^||^ 

< 7r 2 C a 2 ||6||^||/||i a +27T (c a ||a|U^72||fc|U 

< ^Hfcfooll/ll^ + 4vr ((Tr/^qJUall^llftll^ll/ll^ + (0.9062) 2 ||6 
= \\b\H ((0.9062) 2 + 7r 2 C fl 2 + 2 7 r 2 ||a|| 2 C a 4 ) 

using (a + 6) 2 < 2a 2 + 26 2 , the triviality ||/||ia < II/IIh 1 ; and lemma |H1 Therefore 



T" 1 1 "11 OO || J || OO 

2 



2 



hi < 2.2715(1 + ||o|| 0O )C a + 



((0.9062) 2 + vr 2 c7 2 + 2 



7r a 



4x1/2 



□ 



and the claim follows using (a 2 + 6 2 + c 2 ) 1//2 < a + 6 + c for a, 6, c > 0. 

We now turn to the regularity of the eigenfunctions of U. Recall that if Ux = Xx, 
A ^ then formula (fTH|) applies. 

Lemma 10. Suppose k > 1. Suppose a(z),b(z) are analytic functions on a common 
(closed) ellipse E C C with foci ±1 and sum of semiaxes e v = S + s > 1 (see theorem 
HP- £e£ Halloo^ = max^ e s|a(z)| and similarly define ||&||oo,e- Suppose Ux = Xx for 
x G H 1 and A / 0, and suppose \\x\\ni = 1. Let = max^g J\a(e)cZ£ 



A 6(0 



z < 8(sinhr/)" 1 exp (A E + B B /|A|) ke~ kr >. 



max 2gE 
(26) 

Proof. First, formula (|T9*j) determines an analytic continuation x(z) of x(t) defined on I 
to z G E. Apply theorem 4.4 and, specifically, inequality (4.16) in |2Zj to x(z), noting 
that 



max|a;(z)| = \x(— l)|exp ( max 



a(s) + —b(s) ds 
X 



< \x(-l)\exp(A E + B E /\\\) 



and conclude ||x — pWh 1 < \ x (~ l)|8(sinh?y) 1 exp (A E + Be/\X\) ke kv . Also \x(— 1)| < 
0.9062 by lemma El Estimate (JUJ) follows. □ 
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5. Approximating the DDE monodromy operator and its eigenvalues 

5.1. Approximation of U. We propose a Chebyshev spectral collocation method for 
the linear DDE initial value problem (|16j) based upon the ODE method in section El 
Namely, the approximation p G Vn to x G C (I) which solves (fT^j) is defined by 



(27) Dnv = M a v + u, where u 



j 



bit,) fit,), j=0,l,...,A-l, 
/(I), J=N, 



where v G <C N+1 and Dn, M a are (N + 1) x [N + 1) matrices as before. As in lemma 121 
v is equivalent to p: v = EnP, p = Pnv. 

In fact, there is an [N + 1) x [N + 1) matrix approximation to the delayed term 
"b(t)x(t - 2)" in equation (HHJ) defined by 



(M b w) 



b{ti„i)w u 1 <i<N 
w x , i = N + l, 



where w G C N+1 . Compare equation © defining u. 

The collocation method now approximates U itself and not merely solutions to partic- 
ular initial value problems. Let 

(28) U N = (pN-M^Mt, 

a linear operator on C N+1 . (Existence of the inverse in ()28)1 is not automatic, but one 
can show that for continuous a{t) it exists for all N sufficiently large |30j.) 

The finite-dimensional linear map is an approximation to U on H 1 as we will show 
in an a posteriori sense. To be more precise we need to redefine Un to act on H 1 . Recall 
that H 1 C C(I). 

Definition. Let Iljv be the orthogonal projection on H 1 which has range Vn, the degree 
N polynomials. Define G C-iH 1 ) by 

Un = PnUnEn^-n- 

Clearly r&nk(U N ) < N + 1. Noting that r&nge(U N ) C Vn C H 1 , it follows that 
UnP — if and only if p G Vn and Unv = Xv where v = EnP- We will show in theorem 
El that because (in an a posteriori sense) UnJ ~ Uf on a finite-dimensional subspace of 
H , the large eigenvalues of Un approximate those of U. 
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5.2. Eigenvalue perturbation for diagonalizable operators on Hilbert spaces. 

The perturbation result we prove and use, theorem^J generalizes the well-known Bauer- 
Fike theorem for matrices |2], discussed in appendix C. Theorem is probably not new 
and is fairly close to the core idea of the eigenvalue condition number introduced in [3Tj . 
It is also close to proposition 1.15 in [Jj. 

Recall that a separable infinite-dimensional Hilbert space is isometrically isomorphic to 
I 2 (see, e.g. [17 j for this basic fact) where I 2 denotes the space of sequences a = (a 1; a 2 , . . . ), 
a.,- G C, such that |a n | 2 < oo. Denote the standard basis elements of I 2 by 8f (Sj) n = 1 
if n = j and (5j) n = otherwise. Note that A G £(/ 2 ) is diagonal in the standard basis if 
for each j, ASj = Xj5j for some Xj G C. 

Theorem 11. Let X be a separable infinite- dimensional complex Hilbert space (or, in 
the finite- dimensional case, suppose X = C d ). Suppose A G C{X), a bounded operator 
on X , is diagonalizable in the sense that there is a linear, bounded map V : I 2 — > X 
(respectively, V : C d — > X) with bounded inverse and a diagonal (in the standard basis) 
operator A G C{1 2 ) (respectively, A G C dxd ) such that A = VAV' 1 . If B G C(X) and if 
Bx = fix for fi G C and \\x\\ = 1 then 

(29) min \/i - A| < ||^|| HV^ 1 1| || (B - A)x\\. 

xea(A) 

The hypotheses of theorem ^2 imply that o~(A) is the closure of the set of diagonal 
entries of A: o~(A) = Recall that o~(A) C C is compact for any bounded operator 

A G C(X), and thus the "min" in estimate ()29)1 is appropriate. On the other hand, the 
hypotheses do not imply that A is a compact operator. 

Theorem ^2 can be generalized to Banach spaces which are isomorphic to sequence 
spaces (i.e. V spaces). More specifically, we need to hypothesize an operator A, similar 
to A, for which ||(A — /u/) _1 || < sup A . |Aj — /i]" 1 where {Aj} is a dense subset of o-(A). It 
is clear that assuming A similar to a diagonal operator on V suffices, for instance. On 
the other hand, if A is merely normal on a Hilbert space then the just-mentioned 
spectral property will hold. 

Each norm appearing on the right side of estimate ()29)1 is different. Recall that if Vi, 
V2 are normed vector spaces with norms || • ||i, || • H2 and if L : Vi — ► V2 is linear and 
bounded then ||L||v 1 ^v 2 = sup _^ gVl ||Lu|| 2 /||f ||i- We can write the right side of (|2~9"j) in 
a cluttered but precise manner as "H^H^^Hl^ 1 \\x->p \\ (B — A)x\\x" 

Proof of theorem\ll\ If fi G o~(A) there is nothing to prove. Otherwise, (A — fil) 1 = 
V (A — fiiy 1 V^ 1 so that (||(A-/i/) _1 ||) _1 < ||V|| 1 1 1 1 1 / \\(A - fil)' 1 ^ Note (A - 
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fil) 1 is diagonal and bounded with 

II (A — ul)~ 1 \\ < sup | A,- — it I -1 = max I A — Ml 1 = I m i n I A — itl I , 

j \e*(A) \\e<r(A) J 

where {Aj} are the diagonal entries of A. Thus 

ll^llll^ 1 !! 

(30) min |A — u| < - t~t- 

On the other hand, if /i is an eigenvalue of B and x G H 1 satisfies Bx = fix and 
||x|| = 1, then (B — A)x = — (A — fil)x, or x = — (A — fil)~ (B — A)x. Taking norms, 
it follows that 1 < \\{A — fil)~ l \\ \\{B — Combining this with inequality (|3T)|) we get 

estimate (PHJ)- □ 

Recalling lemma El we know that there exist polynomials very close to the eigenfunc- 
tions of U. This motivates the following corollary to theorem ITTl 

Corollary 12. Assume the hypotheses of theorem Additionally, suppose p G X 
satisfies \\x — p\\ < e. Then 



min lit - A| < ||V||||V _1 | 
xea(A) 



\B\\ + \\A\\) + \\(B - A)p\\ 



Proof. Note \\(B - A)a;|| < \\(B - A)(x - p)\\ + \\(B - A)p\\ < (\\B\\ + ||A||)||a; - p\\ + 
\\(B-A)p\\. □ 

5.3. Approximating the eigenvalues of U . Our strategy for estimating the size of 
the difference between the computable eigenvalues of Un (actually, Un) and the desired 
eigenvalues of U (actually, U) follows the outline: 

(i) diagonalize Un = VkV~ l by numerically diagonalizing Un = VAV^ 1 ; 

(ii) estimate \\V\\, ||V^ _1 || by noting PnV = T | . . . | T N V (see lemma US} and 

finding the matrix T numerically; note ||f)jv||ii' 1 < I Ai 1 1| t^|| || V' 1 1| if Ai is the largest 
eigenvalue of Un] 

(iii) also using Un, approximately solve consecutive initial value problems with his- 
tory functions To, T\, . . . , Tn (i.e. the if ^normalized Chebyshev polynomials in 
subsection 14 .2|) and record a posteriori estimates from theorem El 

(iv) use lemma ITU1 to show that an eigenvector of U with eigenvalue away from zero 
is well-approximated by a polynomial p of degree k < N, and thus by a linear 
combination of T , . . . , 5\; 

(v) (this is Theorem 1 14J): use corollary^] with X = H 1 , A = Un, and B — U; estimate 
norm ||?7|| from lemma bound ||(i? — A)p\\ by estimates in (iii); conclude that 
eigenvalues of U are close to the computed eigenvalues of Un- 
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5.4. An a posteriori theorem for the eigenvalues of U . First we show in detail 
how diagonalizing Un diagonalizes Un- 

Lemma 13. Suppose Un = V AV~ l with A = (Aj)^ 1 diagonal. Let t> fc be the kth column 
of'V , so that UNVk = \vk, and definepk = Pn^u- Expandpu in H l -normalized Chebyshev 
series Pk(t) = ^2f=i ^jkTj-i(t), thus defining an invertible matrix T G C( 7V + 1 ) x ( iV + 1 ) . 
Denote by a = (ai, ■ ■ ■ ) an element of l 2 . Let V : l 2 — > H 1 be the operator 

N+l N+l /N+l \ 

fc=l fc>Af+l j=l \k=l / k>N+l 

Then V is bounded and invertible. Furthermore, \\V\\ < (||r||| + l) 1 ^ 2 and \\V~ 1 \\ < 
(l|r _1 ||2 + l) 1 ^ 2 ; where \\ ■ || 2 is the matrix 1-norm on £,( N + 1 ) X ( N + 1 ) . 

Finally, from the diagonal matrix A, define A G C{1 2 ) by (Aa)j = XjOj if j < N+l while 
(Aa)j = otherwise. (Thus A is a diagonal operator of finite rank.) Then Un = VAV^ 1 
and also \\Un\\ < (max|Aj|) || V\\ \\ V~ l \\ . 

Proof. First, V is invertible because V and Pn are invertible and {Tj(t)}jL is a linearly 
independent set in H 1 . Next, the inverse of V is given by 

(oo \ /N+l N+l \ 

j=l J \k=l k=l J 



a map V 1 : H 1 — > I 2 . It is easily checked that V x Va = a for any a G / 2 . 
Next, 



Af+l 



ii^ii^ = E 



AT+l 

E ^fcafc 
fc=l 



2 



/N+l 



+ E N'sK EW + E 

j>7V+l \fc=l / j>iV+l 



< l|r||^||a|| 2 2 + ||a|| 2 2 , 



so IIV^I) 2 < ||r||| + 1, as claimed. A similar calculation shows ||V 1 || 2 < ||T 1 \\ 2 ] + 1. 

Finally, we show Un = VAV^ 1 by action on basis elements of H 1 . If j > N + 1 then 
Ujffj-! = PnUnEnUnTj^ = while VAV' 1 ^ = VASj = 0. For j — 1, . . . ,N + 1, 
however, 

N+l 

(31) C/atT^! = PnUnEnT^ = PnVAV^P^T^ = Y J ^~ 1 )k3 P NVAV~ 1 P N x p k 



k=i 

N+l N+l N+l 

n \ /-n-l\ \ri 

l-l- 

k=l 1=1 ' k=l 



JV+l I\+l 1\+1 

^(T'^kjXkPk = E (E r ^M r_1 )fcj)^- 
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We have used the easily checked facts that ^^/(T ^kjPk — Tj-\ and V 1 P N 1 Pk = &ki 
where is a standard basis element of C N+1 . On the other hand, by definition of V, A, 
and V~\ 

N+1 N+1 

vkv~ l fj_i = v (AiCr-%-, . . . , Aiv+iCr- 1 )^!,,-, o, . . . ) = ^ ( £ r Zfe A fc (r-%)7u, 



i=i fc=i 



exactly the result of (|3T)I . 

Lemma ITB1 defines V" as the product 



□ 



(32) 



V 



T n ... TV T, 



N J- N+1 




where T is an (N + 1) x (N + 1) matrix satisfying 



PnV 



Pi 



Pn+i 



T) 



N 



r. 



Here 



• | . . . | • " denotes a "matrix" whose columns are elements of a function space 
PQ . The right hand factor in equation (1321) is a bounded operator on I 2 with "J" denoting 
the identity operator on the space spanned by {Tj}jL N+1 . In any case, if C is the (N + 
1) x (JV + 1) matrix described in equation (J3J in section El then 

r = v^/2 diag(2" 1/2 , 2, 3, JV + 1)CV. 

We can now state the main theorem on numerical approximation of U. It says that 
the difference between the eigenvalues of the (approximate) monodromy matrix Un and 
the exact eigenvalues of the monodromy operator U can be bounded using the machinery 
developed in this and the last section and also a posteriori information for specific initial 
value problems provided by theorem El 

Theorem 14. Let N > 1. Suppose a,b in DDE (jlfij) are analytic on I = [—1,1] with 
regularity ellipse E D I. Let e v = S + s > 1 be the sum of semiaxes of E. Let = 



max zeE 



f-i> 



B, 



max 2£E 



Suppose fi G C is an eigenvalue of U (and of If) such that |//| > S > 0. Let 

e k = (8(sinhr/)- 1 c Ae+Be/s ) ke~ kri 

for integers k > 1 . 

Assume Un = VAV^ 1 is diagonalizable, with A = (Aj) a diagonal matrix and eigenval- 
ues {Aj}^ 1 ordered by decreasing magnitude. This gives a diagonalization Un = VAV^ 1 
as in lemmaYFh Let cond(V) = \\V\\ ||V^ _1 || • 
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* ~ / k \ 1/2 

Suppose \\UTj — U N Tj\\ H i < Uj for j = 0,...,N and let £ fe — ( ^=0^1) / or ^ = 
l,...,iV. De/me 

(33) u k = e k (\\U\\ + [Ai| condO/)) + (1 + e k )£ k . 
Then 

(34) min \fi — A,| < min {u>i, . . . , u>n} cond(V^). 

j'=l,...,iV+l 

Remark. Note that one computes Vj by solving the ODE initial value problems 

y = a{t)y + b{t)f j {t), y(-l) = f j (l), 

and then by using theoremEl corollary^ and lemma[7| Also, \\U\\ is estimated in lemma 
Elin terms of ||a||oo, \\b\\oc, and C a \ the last constant bounds the fundamental solution 

Proof. Let x G i/ 1 be a normalized eigenvector associated to fi: Ux = fix, \\x\\gi = 1. By 
lemma EH for each k = 1, . . . , N there is q k G 7-^ such that \\x — qkWu 1 < £fe- 

Apply corollary II 21 with ^ = if 1 , A = Un, B = U, and p = q k . Also, use lemma fTTTl to 
diagonalize C/jy = W" 1 , with a(U N ) = {AJf^ 1 U {0} and find 

(35) min \fi-\j\ < 1 1 1> 1 1 I|t>— x || (e k (\\U\\ + \\U N \\) + \\(U - U N )q k \\ w 

j=l,...,N+l \ 



j=l,...,N+l 

On the other hand, q k = YHj=o a kjTj{t) so 



k 



\\(U-U N )q k \\ m < ^2\a kj \\\{U-U N )fj\\ H i < ^2\a kj \ 

3=0 j=0 

< \\qk\\m €k < (IMI/r 1 + hk - x\\ H i)^ k < (1 + 
by Cauchy-Schwarz. From ()35)) and lemma we see that 

min +i \/i - Xjl < cond(y) (e k (||f>|| + |Ai| cond(t>)) + (1 + e fc )&) . 
as desired. □ 

5.5. Using the theorem in an example. The theory of the last section is complicated. 
It remains to show by example that the bounds given there are practical. 

It is worthwhile to give heuristics about the behavior of the constants e k , Uj, £ k , and 
uj k . First, for fixed analytic functions ait) and b(t), e k is monotonic decreasing and 
decays exponentially with rate arbitrarily close to = e~ v < 1. Second, v 3 - is (roughly) 
increasing as j increases. Examples show it is commonly within a few orders of magnitude 
of machine precision e m for j up to approximately 0.7N, supposing N is sufficiently 
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large to achieve e m accuracy for interpolation of the coefficient functions a(t),b(t). The 
meaning of "few orders" is dependent on stiffness, however — see subsection 13.21 Third, 
accumulates the Uj so is also a few orders of magnitude above e m . As N is typically in 
the range 10 < iV < 300, accumulation is small and is at most two orders of magnitude 
above the maximum Uj for j < k. On the other hand, Vj will be large for j pa N because 
UnTj is a poor approximation of UTj in that case, and thus can become 0(1) for 
k pa N. 

We conclude that Uk typically decreases as decreases until (\\U\\ + |Ai| cond(V) j ~ 
£k at which point the minimum value uj is achieved. 

Note, also, that cond(V") is generally an increasing function of N which seems to be of 
size 10 3 to 10 6 for 40 < iV < 300 as we see in the following example. In fact, figure ITTib 
suggests that cond(V) grows subexponentially in the example below. Such behavior is 
evidently desirable for our application. 

Example 5. Consider the DDE, 

(36) y = ay + {b + sin{3irt))y{t-2) 

which is equation (0) of the introduction (section |TJ). Fix parameter values (a, b) = 
(—1.1,1). We find numerically that <j(Un) = = 0.9369 for the monodromy matrix, 
for all iV greater than about 15. We want to confirm this approximation and thereby 
confirm stability for this parameter pair. 

Though b(t) = b + sin(37rt) is entire, we nonetheless must choose a regularity ellipse 
E to find the constants in theorem EI In particular, the oscillation of b(t) for t G J is 
related to the growth of \b(z)\ for z with increasing imaginary part. In fact, note that 

< |6|(1 + S) + (37r)~ 1 (cosh(37rs) + 1). 

Thus one wants to choose the smaller semiaxis s to be relatively small. It turns out that 
the choices s = 0.5 and S = yl + s 2 pa 1.1 are reasonable. 

Figure ITUk shows the growth of the decay of e^, and the (non-monotonic) behavior 
of Uk when = 184. Note that uj = minu;& occurs at k ~ 160. The constants in the 
statement of theorem El are quite large when k is small. In addition to being sensitive 
to the regularity ellipse for b(t), as just described, the lower bound 5 on interesting 
eigenvalues controls the size of iV necessary to achieve a given accuracy. Here we let 
5 = 0.2, that is, we consider the eigenvalues \i of U such that > 0.2. 

The result of applying theorem E] is that the error radius is approximately 0.0434 
so we have proven that < 1 for all eigenvalues of U as shown in figure |2] of section 



B f 



max 



b + sin(37r£) d( 
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Figure 10. (a) Behavior of constants e&, and u)k (see text) in theorem 
d for DDE flU and JV = 184. (b) cond(V r ) as a function of JV. 
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[0 As A" increases above A" = 184, furthermore, we see that the error radius decreases 
exponentially (figure El in section QJ). 



6. Generalization to systems 

6.1. Bounds on fundamental solutions. Although there are other technical compli- 
cations, generalization of the previous results to systems is mainly distinguished by di- 
minished understanding of the fundamental solution to the homogeneous ODE problem. 

Recall, in particular, that theorem El corollary El and lemma El for the scalar case used 
the bound |$ a (£)| < C a = exp (^f\ max{Rea(s), 0} dsj on the fundamental solution 
$ a (t) to y(t) = a(t)y(t). On the other hand, lemma ITU1 required a uniform bound on the 
analytical continuation $\(z),z G E, of $\(t) solving x(t) = (a(t) + b(t) / \)x(t). 

We start with the existence and analytic continuation of the fundamental solution. 

Lemma 15. Suppose A(z) G <C dxd is (entry-wise) analytic on a convex open set E D I = 
[— 1, 1]. Then there is a unique function $4(2), analytic on E, satisfying 

(37) $ A (z) = I d + J* A(0*A«)dC. 

If\A(z)\ < a for z G E then \® A (z)\ < e Q|z+:L| for z G A. Furthermore $ A (t) = A(t)® A (t) 
for tel. 
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Proof. For z G E choose a path 7(5) = — 1 + (z + l)s/r for r = \z + 1| and < s < r. 
Picard iteration on (j3?jl . for z G is 

$,(z) =I d + — f A(7(s))Vi(7(«)) = Id- 

r Jo 

If T]q(z) = \§ q (z) — § q -i(z)\ then induction gives r) q (z) < a q r q /q\, so {§ q (z)} is a Cauchy 
sequence of analytic functions on E. Thus $,4(2) = liniq $ q (z) is analytic. Equation (p!Tj) 
implies = A(t)$ A (t) for i G /. Finally, 

d 



\§a{z)\ = lim 
q 



< ctM/j! = e 

j=0 



(IT 



□ 

One can easily show that, in addition, 

Corollary 16. Fzx s E I. Let Q(t) = for s < t < 1. Then tt satisfies 

Cl(t) = A(t)tt(t), n(s) = I d and furthermore |0(s)| < e 5( *" s) if \A(t)\ < a for r G [s,t]. 

There is, however, a more refined source of bounds on fundamental solutions. Namely, 
one can use the collocation algorithm (below) to approximate the fundamental solution 
and then use a posteriori estimates from theorem El to bound the fundamental solution. 
There is a "bootstrapping" aspect to this: one must have some bound on the fundamental 
solution in order to compute the a posteriori estimates which leads to an improved bound. 
This is the content of lemma below. 

6.2. The collocation algorithm for initial value problems. We now describe the 
collocation algorithm for initial value problems of the form 

(38) y(*)=A(t)y(t) + u(t), y(-l)=y , 

for y(t), u(t), yo G C d and A(t) a d x d matrix. 

As before, fix the interval / = [—1, 1] for t. Define the space of C d -valued continuous 
functions C(J) = C(I) ® C d so that f = (f\, . . . , /^) T is a d x 1 column vector of functions 
fj G C(I). For f G C(J) let ||f ||oo = max tg r \ f{t)\ where "| • |" is the Euclidean norm on 
C d , that is, |v| 2 = Xlj=i \ v j\ 2 - Note | • | induces a norm on d x d matrices; we also denote 
this norm "| • |." Continuous matrix- valued functions A(t) = (aij(t)) by definition satisfy 
G C(J) for all i, j. For such A define the norm \\A\loo = max ie / |^4(£)|- 

Let Vjq be the space of C d -valued polynomials of degree at most N. 

On the function spaces C(7) and Vn we have collocation operators. In particular, 
evaluation at N + 1 Chebyshev collocation points Cn = {to, ■ ■ • , £jv} = {cos(j7r/iV)} to 
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produces a vector in where I = d(N +1): 

(39) S N f = (A(to), • • • , f d (t ), • • • , h(t N ), f d (t N )) T . 

(From now on, we order the scalar components of C' consistently with the output of 
£n-) Restricting Sn to polynomials gives E^v : Vn — > C. The inverse of Ejv, namely 
Pjv : — ► Vn, creates a C d - valued polynomial from collocation values. The interpolation 
operator is I/v = Pn^n '■ C(7) — > Vn- That is, if p = I^f for f G C(I) then p is a C d - 
valued polynomial of degree N such that p(tj) = f(tj). 

The "Chebyshev (spectral) differentiation matrix" is the map 



D 



N 



E N o — o Pjy : 

at 



If Dn is the scalar differentiation matrix (section then D^r = -Djv <8> ^d- In Matlab, 
using cheb.m from j22], if D=cheb(N) then D^r = kron(D,eye(d)). 

Recall that in subsection 13.11 we defined matrices and M a and also a vector u in 
order to write the collocation algorithm as a matrix calculation. Here we define 

\ I u(to) \ 



B N = D N ® I d , M A 



A(t 



N-l, 



\ 



and u 



d J 



u(tiv_i 

\ y / 



where Id, Od are the d x d identity and zero matrices, respectively. Thus Djv, Ma are I x I 
matrices and uGC'. 

Lemma 17 (the collocation algorithm for systems). Let N > 1 and tj = cos(jrj /N), 
j = 0, . . . , N. For A(t) a continuous d x d matrix-valued function of t E I , u(t) G C(I), 
and y G C d , the following are equivalent: 

• p G Vn satisfies 

p(t j ) = A(t j )p(t j ) + u(t j ), 0<]<N-l, and p(-l) = y ; 

• v G C' = C d ^ +1 ) satisfies 

D jyv = M4V + u. 

The equivalence is p = P^v and v = Ejyp. 

Proof. Immediate from the definitions. □ 



For the collocation algorithm we have the following a posteriori estimates: 
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Theorem 18. Let N > 1. Suppose A(t) is a continuous dx d matrix-valued function of 
t e I, u(t) £ C(J), and y G <C d . Suppose 

|$a(*)$a( s ) -1 | <C A for all - 1 < s < t < 1. 

Suppose y G C(J) satisfies the initial value problem (IH8j) . Let p G Pat fre £/ie C d -valued 
polynomial described by lemma\F\ and let R p = p(— 1) — A(— l)yo — u(— 1). T/jen 



(40) 
and 
(41) 



Pi 



< 1C> 



\Ap - IjvCAp)!^ + ||u - /jvHUoo + \TL 



y - p||oo < ^pl^GU + 1) \\Ap - Ijv^p)^ + ||u - /iv(u)||oo + |# F 



Proof. Let q = Jjv(Ap) and w = Jat(u). Since r = p — q — wG Vn and r(t,-) = for 
j — 0, . . . , N — 1, it follows that there is v G C 1 such that r = vl N . (Recall /^(t) is a 
polynomial defined in lemma |3J) Evaluating at t = —1, we find R p = v(—l) N N2 2 ~ N so 
\v\ = \Rp\2 N - 2 /N. 
On the other hand, 



(42) 



p = A(y - p) + (Ap - q) + (u - w) - vZjv, (y - p)(-l) = 0, 



so 



ds. 



y(t) - p(t) = J ^ $a(*)$a00 _1 [(Ap(s) - q(s)) + (u(s) - w(s)) - vfo(s) 
Taking | ■ | norms and using lemma |U to find ||vZjv||oo < , 

|y(t) — p(t)| < 2CU[||Ap-q|| oo + ||u-w|| 0O + |i2 p | 
giving (|4T)|) . For estimate (f4*TJ) . using equation ()42j) and lemmaEl we find that (|4T)|) implies 

SD). □ 

We now return to the computation of bounds on fundamental solutions. For the initial 
value problem (JHHj) . first calculate a typically poor bound |$^(t)$^(s) _1 | < e 2 ^ 00 = 
Ca- Then approximate both the fundamental solution $a(^) and its inverse ^CO -1 by 
the collocation algorithm. Then use the a posteriori bounds on the error in computing 
the fundamental solution to bound the norm of the approximate fundamental solution. 
Repeat as necessary. 

The following lemma gives a precise outline. Recall that ^(t) = $a(£)~ t satisfies the 
"adjoint equation" 4f A (t) = -A(t) T ^ A (t), = Id- 
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Lemma 19. Consider the following initial value problems: 

(43) y s (t) = A(t)y s (t), y s (-l) = e s , s = l,...,d, 

where {e s } is the standard basis for C d . Suppose that for s = 1, . . . , d, ||y 
where p s is the collocation algorithm approximation of y s . 

Note $a(0 — Yi(t) | • • • | Yd(t) is the fundamental solution to y = A{t)y. Let 

®N(t) = [ pi(t) | ... |p d (t) ]■ Ife = Eti^ 2 then \a> A (t)-$ N (t)\ <Uorallte I. 
If, furthermore, ^ A (t) is the fundamental solution to the adjoint equation z = — A(t) T z, 
and if^N^t) is the collocation approximation of^ A (t), and if \^a{P) — ^ jv(^) | < k> for 
t G I, then 

I$a(0^(s) _1 | < (£ + Halloo) (w + ||*jv||oo). 



Proof. 



\$ A (t) - $ N (t)\ = max - $Ar(t)) u| = max 

uGC d ,|u|=l |u|=l 



(y s (t) - p s (t)) 



s=l 



nax > 

u|=l \ ^ 



1/2 



< max 

|u| 



Eiy.(*)-p.(*)i a <e 



by the Cauchy-Schwarz inequality. To conclude, note 
ISaWMs)" 1 ! < |$a(*)||^a(s)| 



< fll$4 -$ 



A — ^N\\oo 1" \\®N\ oo 



AT oo; 



□ 



damped Mathieu equation, corresponds to A(t) 



It follows 



Example 6. The second order ODE x + x + (10 + 9cos(7rt))x = 0, a relatively stiff 

1 
-10 - 9 cos(vrt) -1 

that \A(t)\, \A{t) T \ < \A(t)\ F = ^2 + (10 + 9cos(7rt)) 2 < for all t, where | • \ F 

denotes the Frobenius norm 0. Thus C l = e 2 l' A H- ss 3.5387 x 10 16 is an a priori bound 
on |$a(*)*a(s) _1 |- We use the collocation algorithm with iV = 50 to find ^ N (t),^f N (t) 
approximating $ a (t) , ^f A (t) . The a posteriori estimates from theorem 1181 are computed 
using C A = C x and, as in lemmad we find |$ A (t)$ A (s)- 1 | < C 2 = 2.7117 x 10 9 . This is 
a significant improvement, but also we can now iterate, using C A = Gi in the a posteriori 
estimates to generate C3, and so on. The result is a sequence of bounds 



|$a(0^a(s) _1 | < 3.5387 x 10 16 , 2.7117 x 10 9 , 19.627, 19.587, 19.587, . . . 
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Each number on the right is a bound, with an a priori argument for the first and a posteri- 
ori arguments for the remainder. Evidently they converge superlinearly to about 19.587. 
By looking at (numerical) time- dependent components of the fundamental solutions we 
see that this is near optimal for the Frobenius norm. In any case, this improvement by 
15 orders of magnitude, and comparable improvements for other examples, makes further 
error estimation practical, as we will see in the following subsections. 

6.3. The monodromy operator for a system of DDEs. Let L 2 = L 2 <g> C d be 
the Hilbert space of measurable f = (fi, ■ ■ ■ , fd) T : I ^ C d such that f s G L 2 for 
s = l,...,d. Recall that L 2 = Lj.(I) is the weighted space defined in subsection 14.21 
||/||| 2 = |/(t)[ 2 (l_t2)-i/2 dt Similarly, define H 1 = H 1 ®^, a Sobolev space in which 
first derivatives are controlled. Note that ||f || jj = Ylt=i WfsWjj 1 anc ^ that lemma El says 
|f(^)| — 0.9062 1| f || hi for fixed t G /; our choice of Euclidean norm on C d is important for 
these latter statements. Also, from lemma[7|it follows that ||f|| H i < 2ix ^ 1 1 f 1 1 ^ + ||f||Lj- 
We address the linear DDE 



(44) y(t)=A(t)y(t)+B(t)y(t-2) 

for A, B continuous and periodic (with period T = 2) d x d matrix-valued functions of 
t G I. The monodromy operator for this DDE is U G ^(H 1 ) defined by 



(tf)(t) = <Mt) 



t 



f(l)+ / $ A (s)- L B(s)f (s) ds 



i 



and y n = U n f solves by method-of-steps the initial value problem consisting of ()44j) and 
y(i) = f(t + 2), t G /, f G H 1 . 

If Ux = Ax for A G C then z = Ax solves z = Az + Bx. Thus if A ^ then 
x = (A + A" 1 ^) x. If $ A (t) is the solution to $ A = (A + A" 1 J B)$ A and $ A (-1) = I d then 

X(t) = $A(t)x(-l). 

The next lemmas are systems analogs of lemmas El and fTTH 



Lemma 20. Suppose \$ A (t)$ a(s)~ x \ < C A for all -1 < s < t < 1. Let a 2 = 1 + \\A\ 
and c = 0.9062. Then 

(45) || U || H i < V2~^d (caC A + {{B^ (c 2 + na 2 C 2 A /2) 1/2 ^j . 



2 
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Proof. Suppose f G H 1 and let g(t) = J t _ l ^ A (s)- 1 B(s)f(s)ds. Note ||Uf || H i < \\$ A (t)f(-l)\\ n i+ 
\\&A(t)g(t) || h 1 5 with the first term bounded using lemma 

d 

ll^)f(-l)||^=Ell(^(*) f (- 1 ))*H^ 

k=l 

< 27rX)ll(*A(*)f(-l)) fc ||2o + \UmA(t){(-l))k\\l < 2ndC 2 A a 2 \i(-l)\ 2 . 
fc=i 

On the other hand, | ($ A (t)g(t)) = A(t)® A (t)g(t) + B(t)f(t) so by lemma 

d 

ll<M*)g(f)||iii = 2tE ll(*x(*)g(*))*ll« + (II (A(t)Mt)s(t))k\\oo + ||(5(t)f(t)),||oo) 2 • 
fc=i 

But 

\($ A (t)g(t)) k \ < \* A (t)s(t)\ < I max \<b A {t)<S> A {s)- l \ \B{s)\ \f(s)\ds 

< CaWbw^J 1 \f( s )\ ds < ^cuilsiUlflli,* < yi^HfiiunfiiHi, 

by the Cauchy-Schwarz inequality with weight (1 — s 2 )^ 1 ^ 2 ds. Similarly, 

\(A(t)<f> A (t)g(t)) k \ < yillAHooCAllBlloollfllHi. 

On the other hand, \(B(t)f(t)) k \ < ||-B||oo|f (t)\ < c||-B|U|f || H i. Thus 

\\^A{t)s(t)\\m < 2*d\\B\\l + Z\\A\\lO% + c 2 ) ||f ||^, 

and estimate follows. □ 

Lemma 21. Suppose A, B are analytic dxd matrix-valued functions oft G / with common 
regularity ellipse E C C which has foci ±1 and s^m of semiaxes e v = S + s > 1. Suppose 
Ux = Ax /or A ^ and ||x|| H i = 1. Lei $a(-z) fre £/ie unique analytic continuation of 
$ A (t) /or z & E and suppose C\ is a bound for \$>\(z)\ if z G 7/p = I&x t/ien 

||x - pIU < 8v / rfL7 A (sinhn)~ 1 A;e- fc,? . 

Proof. First, x(z) = $ A (z)x(— 1) is the analytic continuation of x(t), £ G /, to z G -E 1 . 
Furthermore, |x(»| < C A |x(-l)|. It follows from (4.16) of [211 tliat 

d 

ll x " pIIhi < (8C A |x(-l)|(sinhn)- 1 A;e- fcr ') 2 = d (8C A (sinhn)- 1 A;e- fcr ') 2 |x(-l)| 2 . 

3=1 

Because |x(-l)| < 0.9062||x|| H i = 0.9062, we are done. □ 
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6.4. Collocation approximation of the monodromy operator and its eigenval- 
ues. Define 

(B(to) \ 
M B = '• 

B{t N ^) 

where Id, Od are the d x d identity and zero matrices, respectively. The approximate 
monodromy matrix for DDE (jUj) is 



U 



N 



D 



A? 



-1 



M A ) M B eC 



Ixl 



where I = d(N+l), assuming the inverse exists (which we determine by numerical means). 

Let H N : H 1 — > Vn be the orthogonal projection onto Vn, so that if f (t) = Y^JLo &jTj(t) 
then n^vf = ^2f=o a jTj(t). We define the approximate monodromy operator 

tXv = PjvUjvE^IItv G ^(H 1 ). 
The following lemma is the systems generalization of lemma ED in subsection 15.41 

Lemma 22. Let I = d(N + 1). Suppose XJn = VAV -1 with A = (Aj)'- =1 diagonal. Let 
Vfc, k — 1, . . . , I, be the kth column o/ V and define p& = P^v^., a C d -valued polynomial 
of degree N. 

Expand p k in Chebyshev series p k (t) = Y^f=i 12 s =i ^(j-i)d+s,kTj-i(t) ^e s where F i)k G 
C for i, k — 1, . . . , I. The matrix T = (T^) . k=1 is invertible and, in fact, if 

C = (diag(2- 1 / 2 , 2, 3, . . . , N + 1) C) ® I d , 



, U, 2 , . . . J, Uj 



where C is the matrix described by equation then F = CV. 

Let l 2 = / 2 ®C d and denote a typical element by a = (a\, . . . ,af, a\, 
Let V : l 2 -> H 1 5e denned 6y 

JV+l d d 

j=l S=l j>N+l 8=1 

N+l d N+l d d 

= a i r W'- 1 ) <l +» / .(j- 1 ) d + < '^'- 1 (*) e *' + zJ a i^-l(*) e s 

j=l s =l j'=l s'=l i>iV+l s=l 

JTien V is bounded and invertible and furthermore 



|v|| < (l|r|| 2 + i) 1/2 , || v 1 1| < 



v — 1 II 2 



+ 1) 1/2 . 
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Finally, define A G £(1 2 ) by (Aa)y_i)d+ S = Ay_i)d+ S a^ if j < N + 1 and s = 1, . . . ,d. 
If k > I = d(N + 1) i/ien let (Aa)& = 0. (Thus A is a diagonal operator on l 2 o/ ranfc at 
most / .j Then XJ N = VAV -1 and ||Ujv|| < (niaxi<_j<j |Aj|) cond(V). 

Proof. Follow the proof of lemma with appropriate modifications. □ 

We now give our main theorem on the monodromy operator U for DDE (J4"4*J) . 



Theorem 23. Let d > 1 and N > 1. Suppose A,B in DDE (}44j) are analytic d x d 
matrix-valued functions with common regularity ellipse E D I. Let e v = S + s > 1 
be the sum of semiaxes of E. Suppose C\ is an upper bound for \<&\(z)\ for z G E 
and &\{z) the analytic continuation of<&\(t), where <&\(t) is the fundamental solution of 
±=(A(t) + B(t)/\)x(t). 

Let 5 > 0. Suppose Ux = for x G H 1 , ||x|| = 1, and \fj,\ > 5. Let 



(46) e k = 8v / a 7 C , A (sinhA)- 1 A;e" 



kr) 



for k = 1,...,N. Assume XJn = VAV -1 with A = (Aj)' =1 G £(1 2 ) diagonal and 
eigenvalues Aj ordered |Ai| > | As | > . . . . As in lemma, [HH it follows that \Jn = VAV -1 . 
Suppose \\XJ(Tj-i ® e s ) — XJ^iTj-i £g> e s )|| H i < z/| where j = 0, . . . , N and s — 1, . . . , d. 

Let a = EU^tM) 2 and 

uJ k = e k (||U|| + |A X | cond(V)) + (1 + e fc )& 
for k = 1, . . . , N. Then 

(47) min \u — XA < minlcjx, . . . , cjjvi cond(V). 

i=l,„.,d(JV+l) 

Proof. By lemma |^ for each = 1, . . . , N we have ||x — q/clln 1 < e fc where = IfcX&. 
Apply corollary H21 with X = H 1 , A = Ujv, 5 = U, x = x, p = q^, and e = e k . Conclude 

(48) min \fi - A*| < cond(V) (e fc (||U|| + ||U*||) + ||(U - V N )q k \ 

i=l,...,d(N+l) \ 

On the other hand, q fe = ^ =0 S«=i a ljTj(t) ® e s so 



||(u-Ujv)qfe||Hi < 5Z$^KiK < l|q*llHiCfc < (II x IIhi + ||q*-x*llH0& - ( 1 + e fc)^ 

i=0 8=1 



by Cauchy-Schwarz. From (}4*%|) and lemma EH we conclude with estimate □ 
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Remark. As noted is section 16.21 we can produce good a posteriori ( "bootstrapping" ) 
bounds Ca on and |$ y i(t)$ J 4(s)~ 1 | where is the fundamental solution to 

x = A(t)x. This has two consequences for applications of theorem |2H1 First, the a 
posteriori quantities u!j will be small (comparable to machine precision in practice; see 
theorem 118)1 . Second, the estimate for ||U|| will be reasonable (see lemma EST))) . 

On the other hand, the bound C\ on the analytic continuation <&\{t) of the fundamental 
solution §\{t) to x = (A(t) + B(t)/X) x will still be a priori. 4 Such a bound will inevitably 
be big (see the example which follows). Thus tk will be big for small k. Nonetheless the 
essential point of figure El remains: the consistently small values for z/j will allow us to 
choose k < N large enough so that is small. (One must choose iV large enough, of 
course. This can be done using equation (j4T)j) for e k .) 

Said another way, good a posteriori estimates on initial value problems eventually wins 
when competing with large ignorance of the regularity of eigenfunctions of U. 

7. Example: a delayed, damped Mathieu equation 

The DDE x + cx + (1 + cos(7rt))x = bx(t — 2), b G M, c > is a delayed, damped 
Mathieu equation (compare It corresponds to y = A(t)y + B(t)y(t — 2) where 

A(t)J ° M and fl( t )=(° °V 

W \-l - cos(Trf) -cj W \b OJ 

For {b, c) G [—4,4] x [0,3] we have the numerically-produced stability chart of figure HU 
and we expect that the particular parameter pair (6, c) = (1/2, 1) will be stable. In fact, 
for this parameter pair we get the numerical value = 0.5858 for the largest eigenvalue 
of U. The machinery of the previous section allows us to prove the correctness of this 
eigenvalue to two decimal places, in an a posteriori manner, as follows. 

First, using the same method as example El we find an a priori bound |$a(^)| < Ca = 
e 2VE w 134.2 on the fundamental solution $a(£) solving & A = ^(i)^, $a(— 1) — I on 
t G [—1, 1]. We improve this bound by a posteriori iterations as in subsection 16.21 to a 
new bound C' A = 5.12. 

If x is a (C 2 -valued) eigenfunction of U corresponding to eigenvalue A for |A| > 5 > 
then x = (A + B/X)x and \\A + B/XW^ < a/1 + c 2 + (2 + \b\/5) 2 using the Frobenius 
norm. We need to bound the analytic continuation $a(^), for z in a regularity ellipse E, of 
the fundamental solution Q\(t) of x = (A+ B / X)x. Let s = 0.5 and S = vl + s 2 = y/1.25 
be semiaxes of E where E is an ellipse with foci ±1. Then |$a(-2)| is bounded by an a 

4 One might imagine seeking to find &\(z) for z g E by a spectral method. Unfortunately, this problem 
is a system of first order partial differential equations. We don't need an accurate solution to it, in any 
case. Perhaps an a posteriori bound can be sought by some practical method; we don't know of one yet. 
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Figure 11. Stability chart for the delayed, damped Mathieu equation 
x + ex + (1 + cos(nt))x = bx(t — 2). 

priori estimate |$a0)| < exp((l + S)\\A + B/X^) < 4121 = C A for z G E if 5 = 0.3, 
6 = 0.5, and c = 1. 

We now apply theorem |2H1 with (after trying some smaller values) N = 73. The 
constants e^, and cjfc which result are shown in figure IT2*k. The condition number 
estimate for V given in lemma E21 is computed to be approximately 2.898 x 10 8 , and 
this limits our expectations for accuracy of the numerical eigenvalues to perhaps 6 or 7 
digits (given that double precision is about 15 digits). The error radius (the right side 
of equation (|4*7jl and the main result of theorem ESj) is 0.03019. We have the picture 
of eigenvalues in figure 112b. complete with error bounds. The eigenvalues of the exact 
monodromy operator U which exceed 5 = 0.3 in magnitude (shown) are proven to lie 
within the discs of radius 0.0319 (shown) around the numerically computed eigenvalues 
(shown). 

We have proven the stability of this delayed, damped Mathieu equation for the param- 
eter pair (6, c) = (0.5, 1). 

8. Techniques for evaluation of polynomials and norms 

Consider evaluating a polynomial p G Vn described by its values v G C N+1 at N + 1 
Chebyshev collocation points. Concretely, suppose we want to find z = pit) for £ G [—1, 1] 
given Vj = p(tj) at the points tj = cos(jn/N). 

An easy method is described by the single line of Matlab: 

z=polyval(polyf it(tj ,v,N) ,t) 

This method produces considerable error for large iV because polyf it finds the coeffi- 
cients of p in the monomial basis, an ill-conditioned operation on the values v [29\ . 
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Figure 12. (a) Behavior of constants in applying theorem |22] to our de- 
layed, damped Mathieu equation, (b) Two (numerical) eigenvalues of U/v, 
N = 73, of our delayed, damped Mathieu equation exceed 5 = 0.3, as 
shown. Theorem PZH proves that any eigenvalues of the monodromy opera- 
tor U will lie in the discs of radius 0.03019 around the numerical eigenvalues, 
as shown. 

A good alternative is the "barycentric" interpolation of H23J (see also pQ). The following 
simple formula is special to Chebyshev collocation (extreme) points: 



N I N 

p ® = T^T Vj I ^ where 



t-u 



W.; 



-1)72, j = 0oij = N, 
TV, otherwise. 



j=o ' ' j=o 3 

A proof of the numerical stability of this method appears in [T2j; see also jl]. The following 
is a brief Matlab implementation, though not necessarily an optimal one. 
function p=bary (v,t ,N) 

%BARY Barycentric interpolation for Chebyshev points. 

tj=cos(pi*(0:N)/N) ' ; 7, Cheb points 

wj = [1/2; onesCN-1,1); 1/2] . * (-1) . ~ ( (0 : N) ' ) ; 
for k=l:length(t) 

if min(abs (t (k) -t j ) ) > 2*eps % at generic pt 

zz=wj ./(t(k)-tj) ; p(k)=sum(zz.*v)/sum(zz) ; 
else, j=find(abs(t(k)-tj) <= 2*eps) ; p(k)=v(j); end °/ very near Cheb 

end 



Example 7. To estimate the norm ||w — Inu\\oo, one may evaluate \u(t) — p(t)\ at many 
equally-spaced points t in [—1, 1] where p = I^u. Specifically, let u(t) = sin(2t). Figure 
ITBI shows that the evalution of p(t) is done stably by barycentric interpolation up to degree 
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N = 100, but that the "polyval (polyf it ( . . .))" method behaves catastrophically for 
large N. 
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FIGURE 13. Importance of the right polynomial evaluation method: 
Approximations of \\u — InuWoo where u(t) = sin(2t): dots are by 
"polyval (polyf it ( ...))" and circles are by bary.m, above. 



The above is really just a warm-up. The polynomial calculation most used in this paper 
is the estimation of \\p\\oo when p(t) is a polynomial represented by its collocation values, 
or of || o; || oc when a(t) is an analytic function. For estimating such norms there is a much 
better method than evaluating p(t), respectively a(t), at predetermined points. 

Consider the polynomial case. If p(t) = J2k=o a kTk{t) where Xfc(t) = cos(A;arccost) is 
the fcth Chebyshev polynomial then 



(49) 



N 

< ^ \ a k\ 
k=0 



since |Tfe(i)| < 1. Estimate ()49|) requires the coefficients in the Chebyshev expansion 
of p but these are easily calculated by the FFT, as for instance by the Matlab function 

function a = coefft(v) 

"/oCOEFFT Compute Cheb coefficients of an N degree polynomial p(t) 

^represented by its colloc values v_j = p(t_j) where t_j = cos (pi j/N) 
N = length(v)-l; if N==0, a=0; return, end 

v = v(:); °/„ force into column 

U = fft([v; flipud(v(2:N))] )/N; °/„ do t -> theta then FFT 

a = ([.5 ones(l.N-l) . 5] ) ' . *U(1 :N+1) ; 



Thus ||p||oo < sum(abs (coef f t (v) ) ) if v is a vector of the collocation values. 
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If a(t) is analytic on I then we estimate Halloo by approximating a(t) by high degree 
Chebyshev interpolation and then we use (J30J). A very similar issue is addressed in [T]; 
we resolve ours in a similar way, as follows. 

Given a(t), we start with M — 15 (for concreteness) , evaluate at collocation points 
tj = cos(7rj/M), j = 0, . . . , M, and use the FFT to calculate the coefficients of the 
corresponding polynomial. We then inquire if the last few coefficients (say, the last four, 
for concreteness) are small. In particular, if max{|ciM-3|, • • • , \&m\\ < 10e m , where e m is 
machine precision, then we accept the estimate ||o;||oo ~ ||p||oc and use estimate (jjQj) . If 
not, M is doubled (actually, to one less than the next power of 2; this gives efficiency in 
the FFT) and we try again. We stop with M = 2 12 — 1. 

Example 8. Again consider estimating ||w — ijv^Hoo, this time with N = 5 fixed. Thus 
a = u — I5U is the analytic function in question. The following Matlab implements the 
above algorithm: 

b = coefft(u(cos(pi*(0:5)/5))) ; % coeffs of I_N u 

M = 15; last = 8; tol = 10*eps; 
for s = l:last 

tM = cos(pi*(0:M)/M) ; 

a = coefft(u(tM)) - [b; zeros(M-5, 1)] ; % coeffs of u - I_N u 

z = sum(abs (a) ) ; 

if max(abs(a(M-2:M+l))) < tol, break, end 
M = 2*(M+1)-1 

end 

If u(t) = sin(2t), this procedure stops at iV = 31 with an estimate ||w — /5«||oo < 
0.00070975. By contrast, the use of 1000 equally-spaced points in [—1, 1] and bary.m, as 
in 

tp=linspace (-1 , 1 , 1000) ; up=sin(2*tp) ; 

tj=cos(pi*(0:5)/5) ' ; v=sin(2*tj) ; 

for k=l: 1000, INup(k)=bary(v,tp(k) ,5) ; end 

yields ||w — Inu\\oo ~ max(abs(up-INup) ) = 0.00067538 at substantially greater cost 
(roughly 100 times the execution time). 

9. The monodromy operator for the (delay) ^ (period) case 

As noted in the introduction and on page El if the period is equal to (or an integer 
multiple of) the delay, and if the Floquet matrix of a corresponding homogeneous ODE 
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is exactly known, 5 then complex variable methods can be applied to the stability problem 
for linear, periodic, fixed-delay DDEs; see section 8.3 of [H]. In this section we write down 
a variation-of-parameters formula for the monodromy operator in (essentially) the general 
case of one fixed delay and periodic coefficients. We assume no rational relation between 
period and delay. We see that our methods generalize. The same a posteriori results for 
numerical initial value problems and for numerical approximations of the eigenvalues of 
the monodromy operators should be available with sufficient elaboration. 
Consider the system of DDEs (compare subsection Ifi.Hjl : 

(50) y(t) = A(t)y(t) + B(t)y(t-r) 

for y(t) G C d and A, B continuous matrix-valued (i.e. in C dxd ) functions with period T 
exceeding r but less than or equal to 2r: 

0<r<T<2r, A(t + T) = A(t), B(t + T) = B(t). 

Other cases than r < T < 2r will be handled by appropriate modification of the argument 
below. Note T can always be forced to exceed r by choosing the nominal coefficient period 
to be a sufficiently large integer multiple of the minimal period. 

Now suppose f (t) is a history function defined on the interval — r < t < 0. The solution 
of (JSTJj) for one full period T, with initial value y(t) = f(t), t G [— r, 0], is found by solving 
two consecutive ODE problems 

(51) yi = A(t) yi + B(t)f(t-T), yi(0)=f(-r), te[0,r], 
y2 = A(t) y 2 + B(t) y x (t - r) , y 2 (r) = y i (r) , te [r, T] . 

because 2r exceeds T. 

In figureE]we show the schematic map from f to the solution yj_, y 2 on [0, T]. However, 
to be considered as the action of an operator possessing eigenvalues, the solution of the 
initial value problem for system (|5()jl must be redefined to become an operator acting 
from a space to the same space. 

Let y = C[0, T - r] C[T - r, T], that is, let y be C[0, T] broken at t = T - r. We 
denote a typical element of y by (g, f) with g e C[0, T — r] and f G C[T — r, T}. 

The monodromy operator can now be defined as the following composition 

U:(g,f)G^ ig "" g feC[0,T- r] tr Ti atc f G C[-r, 0] 

solve IVP by J^J . r , r evaluate onto [0,T] , , _ A 

' ^ (yi e C[0,r],y 2 G C[r,2r]) i-> (g',?)ey. 

This composition is pictured in figure ITB1 

5 The Floquet matrix for a linear system of non-constant-coefficient ODEs is usually not exactly known. 
Thus the bulk of the paper, and section [S] in particular, addresses a nontrivial case. 
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Figure 14. The initial value problem for (J5Uj) in the case < r < T < 2r. 
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Figure 15. Schematic of action of U on y. 
The "evaluate onto [0,T]" stage perhaps requires explanation. We have 

g'(t)=yi(t), te[0,T-r], and f(t) = 



yi(t), te[T-r,r] 
y 2 (t), te[r,T]. 



Note that though (g, f) need not be in C[0, T] (i.e. g(T — t) need not equal f (T — r)), 
it is still true that ya(r) = yi(r) and thus we may say (g', f) e C[0,T]. In particular, 
the eigenf unctions of U will be continuous functions. 

Furthermore, because the input g is ignored, U has an infinite-dimensional null space 
and thus G cr(U) has infinite multiplicity. Nonetheless, U is evidently compact. Powers 
of U solve initial value problems for equation (jHOjl out to arbitrary t > 0. The nonzero 
eigenvalues of U determine the stability of (f50|) . 
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We propose that a numerical approximation of U can be constructed by using separate 
polynomial approximations of the inputs g and f by interpolation at the (scaled and 
shifted) Chebyshev collocation points in intervals [0, T — r] and [T — r, r], respectively. 
(In fact, we need do no actual interpolation of g.) The ODE initial value problems in (Jo"Tj) 
are then solved by the method of section El so that a posteriori estimates are available 
for the approximations Pi,P2 to yi,y2- Note Pi,P2 are functions on [0, r] and [t,2t], 
respectively. We then interpolate to find qi ~ g' on [0, T — r] and q2 ~ f on [T — r,T]. 
The information in p2 on [T, 2r] is lost by this procedure, a small loss in efficiency if the 
goal is a long-time solution of the initial value problem. (There is no loss of accuracy, 
however, if A, B, f are smooth on (— e, T + e), (— e, T + e), and (— r — e, +e), respectively, 
because yi,y2 connect smoothly at t = r in this case.) The result of these operations 
is a numerical approximation to U acting from the space Vn x © T y N 2 °f polynomials of 
(possibly) different degrees Ni,N 2 back to itself. 

Though we admit that "bookkeeping" complications have kept us from filling out the 
details in these limited pages, the a posteriori and eigenvalue perturbation analysis of 
earlier sections will first give estimates of the accuracy of initial value problems for DDE 
(HIJ) and will furthermore give eigenvalue estimates as well. 

10. Conclusion 

We conclude with a sketch of the "roads not taken" in the current paper and also with 
certain genuinely open problems. 

First, this paper does not address the multiple fixed delays case for linear, periodic 
DDEs. This omission is merely a technical one; all techniques here generalize without 
difficulty to that case. 

Second, we have not yet considered the case of linear, periodic DDEs with nonconstant 
(but state-independent) delays (i.e. x(i) = A(t)x(i) + £>(i)x(t — r(t))), or the case of 
neutral linear, periodic DDEs (M (t)±(t) + Mi(t)i.(t - r) = A(t)x(t) + £(t)x(i - r)), or 
the case of linear, periodic functional equations (x(£) = f^® k(t, s)x(t — s) du(s); [H]). It 
is clear that many techniques in the current paper will generalize if we assume analyticity 
for the relevant functions (t(-), Mj(-), and k(-,s), respectively). We make this statement 
confidently but without precision. 

Now, specific open problems: 

• Under what conditions does U or U (the exact monodromy operator, in any either 
case) actually diagonalize? 

• In fact, do the eigenf unctions of U generically form a Riesz basis ^0] or some 
other (substantially) linearly-independent basis for H 1 ? (Better yet, find an a 
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priori estimate of the condition number of the finite rank operator with columns 
the normalized eigenfunctions of U corresponding to large eigenvalues /i, e.g. such 
that > 5 for fixed 5 > 0.) 
• Given N, estimate ||U — UjvIIh 1 - 

For none of these questions would an answer actually effect the a posteriori methods 
of the current paper. Indeed these questions relate to building an a priori understanding 
of a spectral approximation of monodromy operators for linear DDEs. 
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Appendix A: Constants for theorem El when a(t) 



ao 



To prove the special-case estimates (|1U|) and (jllj) , we exploit cancellation in the second integral 
appearing in We need a lemma which is well-known in the theory of Fourier series, namely, 
that integration-by-parts allows us to estimate the cancellation which occurs in an integral of a 
regular function f(x) against a rapid sinusoid. Actually, we will need the given corollary. 



Lemma. // / G C 1 [a, b] then 



/MfflU 

cos(nx) I 



Proof. Integrate-by-parts. 



< 



n 



(2||/|| 00 + (6-a)||/ / || 00 ). 



□ 



Corollary. Ift*El then 



3 a()( *- s) /jv(s)ds 



-1 



< 



(7r(|a | + 1) +4)e 2 l Rea °l 
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Proof. Substitute s = cos 9 and use lemma 0] to find 



(52) 



o a t fn 



<"""' ' '7.v(-s-! ds 



2 N- 



arccos t 



Let f(9) = e - aocosd (1 - cos(0)) so / > 0, 
lemma 1101 



1 J lloo 



e ~a cose ^ _ cog ^^ sin (_/V0) rffl. 
< 2e IReao| ; and H^li^ < (| aQ | + ^glReaol. By 



f(9)sm(N9)d9 



arccos t 



< l(4|a | + l)+4)el Rea °l 



and the result follows. 



□ 



If Reao < and iV is small then we can sometimes get a slightly better estimate by not 
pursuing cancellation. In fact, if Reao < then 



3 a ^ t - s) l N { S )ds 



< 



< 



\In(s)\ ds 



1 



1 



2 N-1 



2 N-1 

1 — cos 9 d9 



TT 

arccos t 
TT 



(cos(0) - 1) sin(iV0) 



sin(6») 



sin 9 d9 



by lemma 01 

We can now start with equation (|15j) and complete the proof of theorem in this special 
case. Since Jjv(ojj) = aop, and because |$*| = | (t— s) | < m ax{e 2Rea °, 1} for —1 < s < t < 1, 
it follows that 



\y -p\\ < cAu - to I |oo + \(5\ max 



where c\ = max{2e Rea °, 2}. By corollary HU1 and equation (|13|). ||y — p||oo < ci||u — w||oo +C2|i? p | 
for C2 = (7r(|ao| + l) + 4)e 2 l Rea °l/(2iV 2 ). In case Reao > 0, this is the best we can do. Otherwise, 
by the previous paragraph we note ||y — p||oo < 2||n — iu||oo + min{c2, -k / {2N)}\R P \. 



Appendix B: Compactness of the monodromy operator 

In theorem El we use theorem El and an eigenvalue perturbation technique to estimate the 
difference between the eigenvalues of U (eigenvalues of U, actually, but they are the same) and 
the eigenvalues of a matrix approximation to U, as the latter can be found numerically. Of 
course, U is not a matrix and its numerical approximation will (most naturally, anyway) act on 
a space of finite dimension. We do not address, in theorem I14( values in the spectrum of U that 
are not eigenvalues. However, as we show here, U and U are compact and so the only (possible) 
non-eigenvalue spectrum is G C. 

Definition. Let B be a Banach space with norm || • || and denote by C(B) the set of bounded 
operators on B. We say K G C(B) is compact if the image KS of any bounded subset S C B is 
precompact, that is, if every sequence in KS has a convergent subsequence. 

See section 21 of |17| for a proof of the following proposition. 
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Proposition. The set of compact operators is a subspace of C(B) and forms a two-sided ideal 
(that is, if L £ C(B) and K is compact then both KL and LK are compact). Operators with 
finite rank are compact. The spectrum o~(K) of a compact operator K consists of a countable set 
of complex numbers A& £ C. The only possible accumulation point ofa(K) is zero; if dim B = oo 
then £ o~(K). Nonzero elements of o~(K) are eigenvalues. 

Theorem. U £ C(C(I)), defined by equation (|17jl. is compact. 

Proof. Section 22.2 of shows that integral operators (Nkf)(t) = f \ k(t, s)f(s) ds are com- 
pact on C(I) if the kernel k(t, s) is continuous in each variable; this fact is an easy corollary of 
the Arzela-Ascoli theorem. Note that multiplication of / by the bounded functions b, e ±a (' +1 ) 
are each bounded operators; denote by Mb,M±. The operator <5i = / i— > /(l) is finite rank 
and thus compact. Let (Nxf)(t) = f(s) ds, a compact integral operator. It follows that 
U = M + (<5i + NiM-Mb) is compact (because compact operators are a subspace and a two- 
sided ideal). □ 

Corollary. U £ C(H 1 ) is compact. 

Proof. Recall U = U o i H i on H 1 and that the injection i H i : H 1 =— > C(I) is bounded. □ 



Appendix C: The Bauer-Fike theorem for matrices 

Let || • || be a norm for n x n matrices induced by a vector norm on C n . For such a norm we 
define cond(F) = ||V|| ||^ _1 ||, the matrix condition number ofV. 

Theorem (Bauer-Fike [2]). Suppose an n x n matrix A is diagonalized by V , that is, suppose 
A = VAV^ 1 for an invertible nxn matrix V and A = diag(Ai, . . . , A n ), Aj £ C. If B is another 
n x n matrix with eigenvalue \i £ C then 

(53) min - A,| < cond(V)\\B - A\\. 

i=l,...,n 

Proof. Apply theorem fTTI supposing Bx = \xx for x £ C n satisfying ||a;|| = 1. □ 

That is, if A and B are close in norm and A is diagonalizable then any eigenvalue of B will 
be close to one of the eigenvalues of A if additionally V is well-conditioned. More technically, 
the Bauer-Fike theorem says that an upper bound for the condition number with respect to || • || 
of the eigenvalue problem for diagonalizable A is given by the matrix condition number of any 
matrix of eigenvectors V. Note that different diagonalizations of A produce different values for 
cond(y), in general. We suggest the reader explore this theorem in the following exercises. 

Exercise 2. Show that if A is normal then V can be chosen with cond(V) = 1. (Thus the 
eigenvalue problem for normal matrices is well-conditioned.) 
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spectrum of A * 


/ eigenvalue of B 


{radius) <= cond(V) ||B-, 





Figure 16. The Bauer-Fike theorem bounds the distance between eigen- 
values of B and those of A in terms of \\B — A|| and the matrix condition 
number of a matrix of eigenvectors for A. 

Exercise 3. Let A 1 = diag(l,2,3) and A 2 = [1 5 10; 2 -10; 3] (Matlab notation). 
These matrices have the same eigenvalues. Show (numerically), by choosing random but small 
perturbations of the entries, that the eigenvalues of A 2 are much more sensitive to perturbations 
than those of A±. 

Exercise 4. Let A = A 2 in the above exercise. Let AA be a matrix with entries randomly 
chosen from an N(0, 1) distribution. Let B = A + 0.01AA Experiment numerically with 
different diagonalizations V of A and different matrix norms to make (|53[) as close to equality 
as possible. 

Theorem generalizes Bauer-Fike to Hilbert spaces, of course, but also includes the idea 
that B, A do not have to be close in norm to get the result; they need only do close to the same 
thing to x so that \\(B — A)x\\ is small. 
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